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ABSTRACT 


This  thesis  describes  the  development  of  a  computer  program  that 
calculates  the  propagation  loss  for  low  frequencies  in  a  shallow  ocean 
given  the  depth  of  source  and  receiver,  the  sound  speed  profile  of  the 
water,  the  frequency  of  the  source,  and  the  impedance  and  sound  speed  in 
the  bottom.  The  program  does  this,  by  computing  the  sum  of  normal  modes 
for  a  specified  set  of  boundary  conditions.  At  the  surface,  perfect 
pressure  release  is  assumed,  and  the  boundary  condition  at  the  bottom  is 
one  of  impedance  mismatch.  An  effort  was  made  to  develop  a  Fast  Field 
Program,  which  would  use  a  FFT  to  predict  propagation  loss  at  a  variety 
of  ranges  by  solving  for  a  discrete  set  of  wave  numbers,  but  development 
was  not  completed. 
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I.  INTRODUCTION 


A  wide  variety  of  methods  are  available  for  determining  the  pressure 
field  from  a  point  source  in  a  stratified  ocean  of  constant  depth. 
Figure  1.1  shows  that,  by  using  various  transforms,  it  is  possible  to 
derive  expressions  describing  the  pressure  field  in  terms  of  the  linear 
wave  equation,  the  Helmholtz  equation  and,  finally,  the  field  integral. 
Methods  associated  with  these  expresssions  vary,  and  each  has  relative 
advantages  and  disadvantages  in  terms  of  speed,  accuracy,  and  ability  to 
deal  with  complex  situations.  Unfortunately,  no  single  method  is  supe¬ 
rior  in  all  respects  and  there  must  be  trade-offs.  In  some  situations, 
accuracy  might  be  sacrificed  for  speed,  or  the  ability  to  deal  with 
unusual  boundaries  may  override  the  importance  of  speed  or  accuracy. 

Figure  1.2  shows  the  common  methods  of  solving  the  field  integral  to 
determine  the  pressure.  The  Normal  Mode  Solution,  the  Method  of  Stee¬ 
pest  Descent,  The  Multiple  Scattering  Method  and  Direct  Numerical  Integ¬ 
ration  are  "exact"  methods.  The  Fast  Field  Program  (FFP),  a  method  of 
numerical  integration,  must  be  considered  an  approximation,  not  only 
because  it  approximates  a  Hankel  function  by  an  exponential,  but  also 
because  of  its  reliance  on  a  finite  number  of  samples.  Since  Stationary 
Phase  normally  requires  approximations,  then  so  does  its  derivative,  the 
Method  of  Images.  The  Wentzel-Kramers-Brillouin  (WKB)  Integral  normally 
involves  approximations  so  both  Ray  Theory  and  the  WKB  Mode  Equation 
will  seldom  be  exact.  It  is  possible  to  find  exact  solutions  using  the 
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Normal  Mode  Technique,  but  in  this  paper,  as  will  be  seen  in  Chapter  II, 
practical  assumptions  relegate  the  solutions  to  the  approximate  cate¬ 
gory.  This  thesis  is  limited  to  consideration  of  the  Normal  Mode  Solu¬ 
tion,  which  is  computed  by  the  program  EXACT,  and  a  Fast  Field  Program, 
computer  program  FFP. 

Chapter  II  derives  the  field  integral,  then  develops  the  Normal  Mode 
Solution,  and  describes  the  Fast  Field  Program,  FFP.  The  basic  differ¬ 
ence  between  the  methods  involves  the  requirement  of  the  Normal  Mode 
Solution  for  an  accurate  estimation  of  the  horizontal  wave  numbers  for 
each  of  the  modes.  The  FFP  does  not  calculate  the  modal  wave  numbers 
but  computes  the  pressure  contributions  for  a  large  number  of  pre¬ 
defined  horizontal  wave  numbers.  An  expression  to  describe  losses  and 
phase  shift  in  terms  of  the  impedance  and  sound  speed  ratios  between  the 
bottom  and  the  water  is  also  derived  in  Chapter  II. 

Chapter  III  describes  the  programs  EXACT  and  FFP.  EXACT  solves  the 

field  integral  by  first  finding  good  approximations  of  the  mode  wave 

numbers  using  a  matrix  eigenvalue  technique  and  then  refining  them.  FFP 
calculates  the  pressure  field  at  a  discrete  number  of  horizontal  wave 
numbers  without  regard  to  modal  values.  Both  programs  use  the  vector 
sum  of  the  components  of  pressure  from  the  constituent  wave  numbers  to 
calculate  the  propagation  loss  between  the  source  and  the  receiver, 
taking  into  account  absorption  in  the  water,  losses  to  the  bottom,  and 
phase  changes  at  the  surface  and  bottom. 

Chapter  IV  describes  the  testing  procedures  and  the  results.  Within 
this  chapter  the  testing  criteria  are  defined  and  the  adaptive  nature  of 

the  development  process  is  discussed.  Test  results  provided  some  very 
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instructive  lessons  which  gave  an  insight  into  how  normal  modes  can 
describe  pressure  variations  in  both  the  vertical  and  horizontal. 

Chapter  V  discusses  the  results  and  some  of  the  weaknesses  of  the 
progr  am . 
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II.  THEORY 


In  this  discussion  of  normal  modes,  frequency  will  be  very  low  so 
the  absorptive  losses  due  to  the  water  will  be  negligible,  except  at 
long  ranges.  However,  absorption  in  the  bottom  will  often  be  very  high, 
especially  for  the  higher  modes.  The  theory  was  used  to  derive  formulae 
for  the  lossless  case;  then  absorption  was  introduced  as  an  imaginary 
component  of  the  horiztonal  wave  number.  Although  absorption  in  the 
water  is  low  for  low  frequencies,  it  proved  expedient  to  input  abnor¬ 
mally  high  absorption  rates  into  the  test  programs.  By  doing  this,  it 
was  possible  to  verify  that  the  complex  numbers  were  producing  consis¬ 
tent  values. 

The  theoretical  approach  involves  a  series  of  transforms,  some  of 
which  require  approximations  for  simplification.  These  transforms 
successively  take  the  wave  equation  from  the  time  domain  to  the  fre¬ 
quency  domain,  from  three  dimensions  to  two,  from  range  dependent  to 
wave  number  dependent,  and  finally  to  a  function  that  describes  pressure 
as  a  function  of  the  depths  of  the  receiver  and  transmitter,  the  hori¬ 
zontal  wave  number  and  the  horizontal  range.  Normal  mode  theory 
enables  one  to  consider  the  total  pressure  as  the  sum  of  pressures  from 
many  components.  Each  mode  is  stimulated  by  the  transmitter,  or  sensed 
by  the  receiver,  to  a  degree  that  depends  on  the  positions  of  the 
transmitter  and  receiver  in  relation  to  the  modal  depth  function  curve. 
This  curve  is  described  in  Chapter  IV. 
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In  order  to  keep  the  program  simple,  several  assumptions  must  be 
made;  some  limit  the  practicality  of  this  particular  program  and  could 
be  overcome  in  a  more  comprehensive  routine,  while  others  are  of  little 
co nsequence. 

The  main  assumptions,  some  of  which  will  be  discussed  in  detail 
later,  are: 

1.  source  is  cw  and  monochromatic;  it  emits  a  continuous  wave 
at  a  constant  frequency; 

2.  source  is  a  point  radiating  uniformly  in  all  directions; 

3.  medium  is  homogeneous  in  all  respects  in  the  horizontal; 

4.  bottom  and  surface  are  flat  and  parallel; 

5.  surface  is  perfectly  reflecting; 

6.  bottom  can  be  completely  described  by  its  impedance; 

7.  all  energy  transmitted  into  the  bottom  is  lost,  and 

8.  branch  line  integrals  can  be  ignored. 


A.  THE  FIELD  INTEGRAL 

An  omnidirectional  point  source  operating  at  frequency  f(t)  is  posi¬ 
tioned  at  a  depth  z  in  a  water  mass  that  is  uniformly  homogeneous  in 
sound  speed  in.  the  horizontal  plane  and  can  be  described  in  the  vertical 
plane  by  the  sound  speed,  c(z),  at  any  depth. 

As  long  as  the  variations  are  linear,  the  pressure  is  given  by  the 
acoustic  wave  equation: 


VJp  - 


c  2(z)  ‘ 


d^ 

dt5 


-f(t)6(s~s0) 


2.1 


where: 

p  -  p(s,t)  pressure  as  a  function  of  position  and  time, 
f(t)  =*  function  describing  the  source  amplitude, 
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c(z)  =  depth  dependent  sound  speed,  and 
s  =  position  in  Cartesian  coordinates. 

A  Fourier  transform  from  time  to  frequency,  to,  corresponds  to  con¬ 
sidering  the  equation  for  distinct  frequencies  and  yields: 

V2p  +  P  =  -F(to)6(s-s0)  2.2 


where: 

P  =  P(s,to)  a  function  of  position  and  frequency, 

to  =  frequency  in  radians/sec,  and 

F( to)  =  transform  of  the  driving  function,  f(t). 

Because  the  sound  field  is  independent  of  the  azimuthal  angle,  this 
transform  can  be  formulated  in  the  more  manageable  cylindrical  co-ordi¬ 
nate  system.  In  the  new  system  the  wave  equation  appears  as: 


32P  ,  1  3P 

3r 2  r  3r 


3fP 

3z2 


+  k2P 


-F(to) 


<5(z-z0) 

2itr 


dr 


2.3 


where: 

k  =  wave  number,  , 
c(z)’ 

z  =  depth,  and 

P  »  P(r,z,to),  a  function  of  horizontal  range,  depth  and  frequency. 

The  Hankel  transform  allows  one  to  compute  the  changes  of  pressure 
with  a  change  of  depth,  independent  of  the  range.  It  requires  consider¬ 
ation  of  the  vertical  component  of  the  wave  number,  kz: 
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where: 


i!2 

9z2 


+  kzP  = 


_F_ 

2ir 


6(z-z0)  • 


2.4 


p  =  p(£,rw)  a  function  of  horizontal  wave  number,  horizontal  range  and 
frequency, 

kz  =  k2  -  £2 

£  =  horizontal  wave  number,  and 

k  =  wave  number,  f  .  . 

c(z) 

This  equation  is  then  solved  by  variation  of  parameters.  Let 


p  =  w^zJuCz)  +  w2(z)v(z)  2.5 


where: 


/  \  F 

«,(z)  .  -yj 


fZ 


v(z)<5(z-z0)dz 


W, 


2.5a 


uv 


W,(z)  =  ■7T— 


F_ 

2  IT 


u(z)<5(z-z0)dz 


W. 


2.5b 


uv 


Wuv  =  u  ^  -  v  ^  (the  Wronskian),  2.5c 

and  both  u  and  v  are  solutions  of  the  homogeneous  form  of  equation  2.4. 

Mow  if  v(z )  satisfies  the  upper  boundary  condition  and  u(z)  satisfies 
the  lower  boundary  condition,  then: 
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w 


1 


-F  v(zq)  z  >  z o 

,2*  wuv 

^0  z  <  z0 


2.6a 


and 


w2 


’  F  u(zo) 
,  2  IT  Wyy 
0 


z  <  z0 
Z  >  z0 


2.6b 


This  means  that  the  transformed  pressure  can  be  written  in  terms  of  depth 
dependent  u(z)  and  v(z)  as: 


P  = 


_F_ 

2tt 


u(z0)v(z) 

Wuv 


for  z  <  z0 


2.7a 


P 


p  u(z)v(z0) 

2*  VW 


for  z  >  z0 


2.7b 


A  convenient  way  of  representing  this  is 


u(z>)v(z<) 


2.8 


where: 

z>  represents  the  larger  of  z  and  z0,  and 
z<  represents  the  smaller  of  z  and  z0. 

To  obtain  the  actual  pressure  one  must  next  take  the  inverse  Hankel 
transform  of  zero  order: 
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p  = 


pJ0(£r)CdC 


2.9 


F_ 

2tt 


u(z>)v(z<) 


W. 


JoUDCdC 


2.10 


uv 


Use  is  made  of  the  fact  that  [Ref.  1  ] 


J0Ur)  =  1/2 


Ho^tSr)  +  Hq2)(^) 


2.1  1 


and 


Hj2)(Cr)e-^  =  -Hq1  ^(  £r ) 


2.12 


It  is  then  possible  to  express  the  field  integral  (Equation  2.10)  in 
terms  of  a  Hankel  function  of  the  first  kind  [Ref.  2]: 


P 


oo 


iF 

4-tt 

J  “  oo 


u(z>)v(z<) 


H^1)(Cr)5d5 


2.13 


The  Hankel  function  can  be  approximated  by 


(CD 


l/-§-  i(5r-n/i0 

V  iter 


+ 


_j _ 

18  £r 


+■ 


2.14 
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Restricting  this  .expansion  to  the  first  term  makes  the  integral  much 
more  manageable  but  restricts  the  minimum  horizontal  range  so  that  £r 
>>  1  : 


P 


F 

7T 


r  oo 


-00 


u(z>)v(z<) 


e15r^d5 


2.15 


This  field  integral  represents  the  pressure,  within  the  limits 
imposed  by  truncating  the  expansion  of  the  Hankel  function.  This  integ¬ 
ral  forms  the  basis  for  the  derivation  of  formulae  describing  pressure 
using  two  different  numerical  approaches.  At  this  point,  the  approach  to 
the  normal  mode  solution,  which  will  be  used  in  the  program  EXACT,  takes 
a  different  tack  from  the  approach  which  will  be  used  in  the  FFP  solu¬ 
tion  and  program. 


3.  NORMAL  MODE  METHOD  FOR  PRESSURE  DETERMINATION 

In  this  paper  the  normal  mode  method  refers  to  the  process  whereby 
approximations  for  the  eigenvalue  of  a  mode  are  successively  refined  by 
numerically  integrating  the  square  of  a  depth  function  across  the  depth. 
A  correction  term  is  developed  from  this  integral  and  is  used  to  refine 
the  approximation  until  the  correction  is  negligible. 

Let  the  guess  at  the  normal  mode  horizontal  wave  number,  £n,  be  £. 
For  5,  v(z)  satisfies 


v"  +  k|  v  =  0 


2.16a 
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where: 


k 


2  '  _ 
Z  - 


k2  -c,2 


For  5n,  vn(z)  satisfies 


*zn 


v  =  0 


2. 1  6b 


where 


kin  =  k2  -5n2- 

Multiplying  2.1  6a  by  vn  and  2.1  6b  by  v  and  subtracting  yields 


vnv"  “  vvn"  +  vvn(kl"k|n)  =  0 


2.17 


Then 


_d_ 

dz 


(vnV  ~vvn' ) 


+  vvnCk2-^2-^^2) 


0,  or 


_d_ 

dz 


(vnv* -vvn' )  =  U2"^)  vvn* 


2.18 


Integrating  from  depth  0  to  depth  h: 


(vnv' 


U2  - 


vvn  dz 


2.19 


or 
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vn(h)v'(h)  -  vn' (h)v(h)  -  vn(0)v'(0)  +  vn(0)v(0)  =  (52_5n) 


h 


vvndz.  2.20 
JO 

The  following  argument  holds  for  any  set  of  permissible  boundary 
conditions  but,  in  order  to  clarify  the  discussion,  the  specific  case  of 
pressure  release  at  z  =  0  and  rigid  at  z  =  h  will  be  selected.  The 
function  v(z)  is  not  a  mode,  but  it  can  be  forced  to  satisfy  one  of  the 
boundary  conditions  without  losing  generality.  If  v(z)  ia  defined  to 
satisfy  the  boundary  condition  at  z  =  0,  then,  for  any  permissible  boun¬ 
dary  condition,  the  last  pair  of  terms  on  the  left  hand  side  of  equation 
2.20  is  identically  zero.  Therefore, 


vn(h)v’(h)  -  vn' (h)v(h)  =  U2  -5*) 


vvn  dz. 


2.21 


For  the  simplest  case  considered  in  this  paper,  the  bottom  boundary  con¬ 
dition  is  rigid  so  that  vn’(h)  =  0.  Therefore 


vn(h)v’(h)  =  U2 


rh 


vvn  dz 

0 


2.22 


and,  since  vz  vn  for  E,  close  to  £n, 


AS2  “C2"^ 


v(h)v'(h) 

/v2dz 
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5n  =  S2  “ 


v{h)vT (h) 
/v2  dz 


2.23 


Besides  proving  useful  in  determining  successively  better  approximations 
of  5n,  this  same  integral  forms  part  of  the  expression  for  the  contribu¬ 
tion  of  a  mode,  since  it  is  proportional  to  the  derivative  of  the  Wron- 
skian  with  respect  to  wave  number. 

The  Cauchy  integral  enables  one  to  write  equation  2.13  as  the  sum  of 
the  residues  of  all  the  poles  of  the  integrand.  Each  pole  corresponds 
to  a  normal  mode. 


P 


iF_ 

2  7T 


7 


unvn^0  ^nrHn 


dW 

dS 


Cn 


2.24 


where: 


WUn)  =  0 

To  derive  an  expression  for  (dW/d£)  multiply  both  sides  by 

un(h)  un' 

B  =  Vh)  (6  3130  equals  ^7  } 


un(h)v* (h)  -  u/1(h)v(h)  =  (  C~Cn)U  +  Cn)  B 


v2dz 


2.25 


The  left-hand  side  of  this  equation  gives  the  Wronskian  for  u  and  v, 
and,  as  ^  ■*  £n,  W(£n)  -»  0,  as  it  should: 
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dW 


2.26 


W(0"W(5n) 

5-Cn 


6UHn) 


•h 

Jo 


v2dz 


2.27 


Since  £  _  £n,  and  by  properly  constructing  the  program,  the  constant  8 
can  be  forced  to  be  equal  to  1, 


dW 

d£ 


Sn 


rh 


2Sn 


v2dz 

•  0 


2.28 


and  equation  2.24  can  now  be  written: 


P 


n 


unvn 


2^n/v2dz 


^0  ^5nr)5n- 


2.29 


Using  the  exponential  for  the  Hankel  function,  as  in  eqn.  2.13: 


iL  y  J_2_  unvn  iCnr  -i  tt/4 
4tt  ^  y  ir£nr  /v2 dz 


2.30 


n 


un  vne 


Hnr 


i  ti/4 


/v2  dz 


2.31 
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In  reality,  this  only  represents  the  pressure  in  a  system  for  which 
there  are  no  boundary  losses:  i.e.,  for  self-adjoint  boundary  conditions. 
If  real  losses,  including  scattering,  were  to  be  taken  into  account, 


fh 


the  term 


JO 


v2dz  must  be  replaced  by  [Ref.  3]: 


•h 

■  0 


v2dz 


v2(0) ( T?)  "  v2(h)(^) 
df;  d£  _ 


2.32 


where: 


.  dv/dz  du/dz 

a  =  -  y  =  - 

v  u 

Ignoring  boundary  losses,  equation  2. 31  represents  the  total  pressure 
at  a  point  which  is  at  depth  z  and  horizontal  distance  r  from  a  point  cw 
source  at  depth  z0.  This  equation  is  the  basis  of  the  computer  program 
EXACT,  which  is  discussed  in  detail  in  Chapter  III. 

C.  DIRECT  EVALUATION  OF  THE  FIELD  INTEGRAL 

The  second  method  of  determining  pressure  at  a  point  requires  calcu¬ 
lating  the  field  integral  itself  and  then  numerically  integrating.  Equa¬ 
tion  2.15  represents  the  pressure  field  in  terms  of  the  Wronskian  and 
the  source  and  receiver  depth  functions  u(z)  and  v(z).  The  integral  is 
solved  by  replacing  it  with  a  summation  across  an  interval  of  wave 
numbers.  Pressure  at  a  point  is  found  to  be  the  sum  of  the  pressure 
contributions  from  each  horizontal  wave  number: 
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^mAC  A? 


2.33 


F 

4tt 


“  ITT  /  4 
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5n 


-Sn 


u(z>)v(z<) 


gimACr 


u(z>),  v(z<),  and  Wuv  are  calculated  from  the  known  boundary  conditions 
at  the  surface  and  at  the  bottom.  AC  must  be  small  enough  to  ensure 
accuracy. 

It  is  impractical  to  pursue  this  theme  without  introducing  absorp¬ 
tion.  Evaluation  of  equation  2.31*  when  mAC  =  Cn  produces  Wuv  =  0  in  the 
denominator,  leading  to  an  infinite  solution.  The  magnitude  of  Pn  then 
depends  on  how  near  mAC  comes  to  Cn-  Introduction  of  absorption  as  an 
imaginary  component  of  horizontal  wave  number  displaces  the  poles  from 
the  real  axis  and  ensures  that  the  Wronskian  does  not  go  to  zero  along 
the  path  of  integration.  This  has  the  effect  of  dispersing  the  energy 
out  of  the  theoretically  limitless  pressure  function  that  results  when 
the  Wronskian  goes  to  zero  at  the  exact  wave  number  corresponding  to  the 
normal  modes. 

The  terms  within  the  summation  sign  of  equation  2.33  have  the  same 
form  as  the  discrete  Fourier  transform,  normally  considered  for  time,  t, 
and  frequency,  oj: 

N-1 

y~  G(nAoj)  gi^unm/N  _  2.34 

n=0 


By  simply  renaming  the  variables,  the  transform  becomes: 
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N 2}  .  ‘ 

G(nA 5)  el2linra/N  =  g(mAr).  2.35 

n=0 

This  is  the  basis  for  the  technique  known  as  the  Fast  Field  Program 
[Ref.^l].  By  comparing  equation  2.3^  to  equation  2.35,  the  following  can 
be  seen  to  be  true: 


G(nAO 


u v 
W 


The  pressure  at  range  R  can  then  be  written  as 


2.36 


PR  “  ^mAr 


_F_ 

4tt 


irmAr 


-iir/4 

e 


A5CFFT[G]]. 


2.37 


The  capability  of  this  Fast  Field  Program  will  be  limited  by  the  number 
of  range  (and  wave  number)  increments.  The  wave  number  sampling  incre¬ 
ment,  A5,  the  range  resolution,  A r,  and  the  number  of  points  in  the  FFT 
can  be  seen  to  be  related  by: 


2irnm 

N 


£r  =  (nA£)(mAr) 


2.38 


or 


2tt 

N 


A£Ar, 


2.39 
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which  results  from  comparing  equations  2.3^  and  2.35.  The  number  of 
sampling  points,  N,  places  practical  limits  on  either  the  maximum  range 
or  the  range  resolution. 


D.  ABSORPTION  AND  BOUNDARY  EFFECTS 


Sound  energy  is  "lost"  by  absorption  in  the  water  and  transmission 
into  the  bottom.  In  both  cases  a  primary  variable  is  frequency  but  the 
angle  of  travel  of  the  mode  and  the  type  of  bottom  must  be  taken  into 
account.  This  paper  limits  discussion  to  a  non-scattering  surface,  water 
volume  and  bottom.  It  also  neglects  energy  that  is  transmitted  into  the 
bottom  and  then  returned  to  the  water  after  refraction  within  the 
bottom. 

Sound  absorption  in  the  water  at  low  frequencies  (under  1  KHz)  is 
due  primarily  to  relaxation  of  the  boric  ion.  Standard  texts  [Refs.  5,6] 
define  absorption  in  this  frequency  range  as: 


0.1F2 

1  +F2 


2v  40 


where: 

a  =  absorptive  loss  in  dB/m,  and 
F  =  frequency  in  kilohertz. 

This  loss  can  be  converted  to  the  fractional  loss,  a,  for  a  given  dis¬ 
tance  in  nepers/meter: 


Loss  =  20  log-jo  eCtr  dB 
ar  =  20  log-]Q  ear 
=  20ar  log-]0  e 
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=  8.7  or 


a  =  8.7a  dB/meter,  or 


=  g— -  nepers/meter 


2.42 


This  loss  rate  can  be  introduced  into  the  equations  derived  earlier 
by  adding  absorption  as  an  imaginary  component  of  the  wave  number: 


kc  =  k  +  ia 


2.42 


The  exponential  form  provides  the  absorptive  loss  factor: 

ikcR  i(k+ia)R 
e  u  =  e 


ikR 

=  e 


~aR 


e 


2.43 


However,  all  of  the  expressions  for  pressure  are  given  in  terms  of  the 

£  r 

horizontal  range,  r,  not  actual  "path"  length,  R.  Since  —  =  the  loss 
factor  can  be  written  as: 


loss  =  e~a  E,  r.  2.44 

Bottom  interactions  cause  losses  and  phase  changes  which  must  be 
taken  into  account.  In  their  discussion  of  propagation  loss  using  normal 
modes  ands  an  impedance  boundary  condition,  Koch  and  Lindberg  [Ref.  7] 
begin  their  discussion  by  defining  the  relationship  between  pressure  and 
its  derivative  with  respect  to  depth: 
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dP 

dz 


2.45 


ikz 


1-R 
1  +R 


P, 


where 

P  =  fluid  pressure  in  the  water  at  the  boundary, 
dP 

—  =  derivative  of  pressure  with  respect  to  depth, 
kz  =  vertical  wave  number,  and 
R  =  complex  relfection  (Rayleigh)  coefficient. 

The  Rayleigh  reflection  coefficient  for  the  boundary  between  two 
homogeneous  fluids  is  given  [Ref.  5,  Equation  6.30]  by: 

Pbcb  _  s^eb 

a  -  2.«6 

Pbcb  +  sineb 

Pw°w  sin0w 

where 

cb  =  sound  speed  in  the  bottom  at  the  interface, 

cw  =  sound  speed  in  the  water  at  the  interface, 

pw  =  density  of  water, 

pb  =  density  of  the  bottom, 

0W  =  angle,  measured  from  the  horizontal  at  which  the  wave  strikes 
the 

bottom,  and 

0b  =  angle,  measured  from  the  horizontal,  at  which  the  wave  is 
transmitted  into  the  bottom. 

Substituting  this  value  of  R  into  Equation  2.45  and  simplifying: 
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dP 

dz 


P. 


-  ikz 


P  wGw 
PbGb 


sin0b 

sin0w 


2.47 


But 


sin©)-,  =  / 1  -  cos 2  05 


2.48 


Snell's  Law  relates  the  incident  and  transmitted  angles  by: 


cb 

cos  ©b  =  —  cos  0W 
cw 


So, 


sin  0b  * 


cos20w  . 


2.49 


2.50 


The  critical  angle,  0C,  is  defined  as  that  grazing  angle  for  which 
the  transmitted  angle  is  zero.  This  means  that: 


G0S  8w  _  cos(O) 
Gw  Gb 


and 


Gb  =  1 

cw  cos0c’ 


2.51 


where  0C  is  the  critical  grazing  angle,  measured  from  the  horizontal. 
Therefore,  from  Equation  2.50, 


sin0b  = 


cos20w 

COS  2  0C 


2.52 
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By  definition, 


cos  ew  =  £ 

and 


kz 

S1^  0w  =  ~  • 


Substituting  Equations  2.52  and  2.54  into  Equation  2.47: 


2.53 


2.54 


dP  _  .  P wcw  /  _  cos26w 
dz  ~  1  PbCb  y  cos20c 


2.55 


In  terms  of  wave  numbers,  Equation  2.53  can  be  used  to  give: 


or 


dP 

dz 


dP 

dz 


P  wc 
Pbcb 


i 


P  wcw 
Pbcb 


2.56 


2.57 


Equation  2.56  shows  that,  for  any  grazing  angle  less  than  critical, 
the  pressure  derivative  will  be  real  whereas,  for  grazing  angles  greater 
than  critical,  the  derivative  will  be  purely  imaginary.  If  the  pressure 
were  complex,  the  derivative  would  be  complex  and  complicated,  too 
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complicated  to  solve  in  all  cases  with  the  basic  program  design.  Therefore, 
only  the  real  component  of  pressure  was  used  in  the  determination  of  the 
boundary  conditions. 

A  practical  discussion  of  the  physical  significance  of  the  reflection 
coefficient  and  the  relationship  amongst  it,  the  grazing  angle  and  criti¬ 
cal  angle  is  included  in  Chapter  IV. 
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III.  PROGRAMS 


Implementation  of  theory  into  the  computer  programs  required  two 
distinctly  different  methods.  The  program  EXACT  made  accurate  estimates 
of  the  normal  mode  horizontal  wave  numbers  in  order  to  compute  the 
pressure  as  the  sum  of  the  pressure  components  from  each  of  the  modes. 
Absolute  values  of  pressure  are  of  no  concern  in  either  EXACT  or  EFP 
because  the  purpose  is  to  compute  propagation  loss.  It  was  assumed  that 
forcing  function,  F,  had  a  value  of  unity.  The  calculated  pressure  can 
then  be  compared  to  the  pressure  at  one  meter,  and  from  that  ratio,  the 
propagation  loss  is  calculable  using  decibel  ratios.  The  program  FFP  did 
not  compute  horizontal  wave  numbers  of  modes;  it  assumed  a  large  number 
of  wave  numbers  equally  spaced,  calculated  the  pressure  for  each,  and 
summed  the  individual  contributions  to  give  a  final  pressure. 

Both  programs  were  written  using  complex  numbers  and  double  preci¬ 
sion  (64  bit  words),  except  for  subroutine  'eigens*  which  was  restricted 
to  single  precision  because  it  relied  upon  IMSL  routine  EQRTIS  which  was 
available  only  in  single  precision. 

Subroutines  f  modes  1f  and  ' modes2f  provide  many  common  functions  for 
both  programs  EXACT  and  FFP.  They  calculate  the  transmitter  and 
receiver  modal  pressure  functions  by  stepping  upward  or  downward  from 
one  of  the  boundaries,  computing  the  pressure  and  its  derivative  with 
respect  to  depth  by  means  of  a  fourth  order  Runge-Kutta  technique 
[Ref. 7].' 
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A.  NORMAL  MODE  PROGRAM 


This  program  is  attached  as  Annex  B.  EXACT  computes  the  approximate 
modal  horizontal  wave  numbers  by  means  of  the  subroutine  'eigens'.  It 
then  refines  the  horizontal  wave  number  estimate  to  the  desired  accuracy 
using  subroutine  'modes!'  or  'modes2'  and  computes  the  pressure  factors 
for  each  mode.  Its  aim  is  to  compute  the  total  pressure  and,  hence,  the 
propagation  loss  at  each  range  of  interest. 

Subroutine  'eigens'  uses  a  matrix  method  [Ref. 8]  to  approximate 
all  the  real  modal  horizontal  wave  numbers  as  well  as  the  first  evanes¬ 
cent  mode.  At  the  ranges  of  concern,  it  is  unlikely  that  more  than  the 
first  evanescent  mode  would  be  significant.  The  matrix  of  wave  number 
estimates  is  found  in  a  relatively  short  time  and,  although  not  guar¬ 
anteed  to  be  accurate,  it  is  complete;  no  modes  are  skipped. 

The  next  step  involves  finding  the  exact  values  of  the  horizontal 
wave  number  for  each  mode.  The  procedure  utilizes  the  fact  that,  for 
each  mode,  both  boundary  conditions  for  a  depth  function  (u  or  v  from 
Chapter  II)  will  be  met.  By  setting  the  depth  function  or  its  derivative 
with  respect  to  depth  to  a  known  boundary  condition  at  the  top  or  the 
bottom,  it  is  possible  to  find  a  vertical  wave  number  that  will  result 
in  the  boundary  conditions  being  met  at  the  the  other  boundary.  The 
boundary  with  higher  sound  velocity  will  experience  an  exponential  decay 
in  the  depth  function.  By  starting  at  this  boundary,  it  is  a  relatively 
simple  matter  to  solve  for  vertical  wave  numbers  that  result  in 
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sinusoidal  variations  that  lead  to  solvable  boundary  conditions  as  the 
other  boundary  is  approached.  This  avoids  the  possible  problem  that 
would  exist  if  an  attempt  had  been  made  to  start  at  the  boundary  which 
had  the  lower  sound  velocity.  In  that  case  it  is  possible  that  instead 
of  solving  for  an  exponentially  decaying  function,  the  function  may  be 
exponentially  growing,  and  no  solution  is  possible.  Figure  3.1  illus¬ 
trates  how  a  poor  choice  of  starting  boundary  could  lead  to  an  unsolvable 
situation.  For  each  of  the  modes  in  the  'main*  routine  calls  either 
’modest1  (when  sound  speed  is  greater  at  the  top  than  at  the  bottom) ,  or 

’modes2’  otherwise.  It  calls  with  starting  estimates  of  the  modal 

horizontal  wave  number;  the  value  from  ’eigens’  the  first  time  and  an 

improved  value  each  subsequent  time. 

*  Modesl  f,  using  the  best  estimate  for  £n,  computes  the  modal  depth 

function  at  each  increment  of  depth,  starting  at  the  surface  where  the 

depth  function  is  zero  and  its  derivative  with  respect  to  depth  is  set  to 

1.0.  The  integral  (Equation  2.31)  is  calculated  at  each  interval  and  its 
final  value  at  the  bottom  is  used  in  the  correction  term. 

At  the  bottom,  the  depth  function  and  derivative  are  compared  to  the 
boundary  conditions  expected  for  a  mode.  For  a  rigid  bottom,  a  deriva¬ 
tive  of  zero  is  expected.  For  an  impedance  bottom,  the  expected  deriva¬ 
tive  will  be  defined  by  Equation  2.56.  The  error  in  the  derivative,  tog¬ 
ether  with  the  depth  function  and  the  integral  is  used  to  calculate  a 
correction  for  the  modal  horizontal  wave  number. 

If  the  surface  sound  speed  is  greater  at  the  bottom  than  at  the  top, 
subroutine  ’modes2’  is  called.  It  performs  the  same  function  as  'modesl1 
but  starts  at  the  bottom,  where  the  depth  function  and  derivative  are 
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( b  J 


Figure  3.1.  Possible  Solutions  to  Depth  Function  Starting  at: 

a:  Boundary  with  Exponential  (Convergent)  and 
b:  Boundary  with  Sinusoidal  Changes  (Possibly  Divergent). 


36 


set  to  meet  the  boundary  conditions  for  an  impedance  bottom  (Equation 
2.56).  A  correction  is  made  after  the  subroutine  has  stepped  to  the 
surface  and  comparison  has  been  made  with  modal  boundary  condition 

expectations;  in  this  case  depth  pressure  function  equal  to  zero. 

'Modes'!'  or  ’modes2'  will  be  repeatedly  called  until  £n  has  been 

satisfactorily  estimated.  The  'main'  routine  then  computes  the  'incom¬ 
plete'  modal  pressure  from  the  modal  depth  functions,  the  integral,  and 
the  horizontal  wave  number,  £n,  (Equation  2.31)  for  each  mode,  storing 
them  for  later  use.  It  is  termed  'incomplete'  because  there  is  no  range 
dependency  at  this  point. 

After  the  last  mode  has  been  processed  by  '  modes  1'  or  ’modes2', 
range  is  introduced  and  partial  pressures  for  each  mode  are  calculated, 
keeping  an  updated  sum  (complex)  of  all  the  modal  pressure  factors. 

These  pressure  factors  are  then  converted  to  propagation  loss.  It  is  a 
simple  matter,  once  all  the  eigenvalues  have  been  found,  to  compute  the 
propagation  loss  for  a  series  of  ranges,  as  can  be  seen  in  Appendix 

B,  Program  EXACT. 

B.  FAST  FIELD  PROGRAM 

This  program  was  written  with  the  intent  of  completing  a  Fast 
Fourier  Program,  FFP,  which  would  utilize  a  Fast  Fourier  Transform,  FFT, 
subroutine.  However,  the  program  was  only  completed  to  the  point  where 
pressure  factors  were  computed  for  individual  ranges  (as  in  EXACT), 
rather  for  a  large  number  of  ranges  (as  is  the  intent  of  an  FFP).  The 
computer  program  FFP,  although  incomplete,  is  included  as  Appendix  C. 
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It  was  mentioned  previously  that  EXACT  and  EFP  are  provided  many 
common  functions  by  rmodes1f  and  1 modes2\  There  are  also  some  basic 


differences 

in  the  way  the 

subroutines  are  used 

by  the  two 

programs. 

In 

FFT,  since 

no  exact  wave 

numbers 

ar e  to  be 

calculated, 

there  is 

no 

requirement 

for  subroutine 

r  eigensT. 

For  the 

same  reason 

,  there  is 

no 

requirement 

to  find  corrections  to 

the  wave 

numbers,  so 

' modes  1 ' 

and 

fmodes2T  are  used  differently.  In  fact,  both  subroutines  are  required  by 
FFP. 

Using  the  incremental  wave  number,  nA£,  subroutine  !modes1!  starts 
at  the  surface  with  the  same  boundary  conditions  as  in  EXACT  and  steps 
down,  computing  the  pressure  function  and  the  derivative,  through  the 
depth  of  the  upper  of  the  transmitter  or  the  receiver,  and  stops  at  the 
lower  of  the  two.  Values  are  required  at  the  lower  level  for  later  use 
with  the  values  from  Tmodes2T  to  calculate  the  Wronskian.  Then  fmodes2! 
starts  at  the  bottom  with  the  same  boundary  conditions  as  EXACT  and 
steps  upward  until  the  pressure  function  and  derivative  have  been  calcu¬ 
lated  for  the  lower  level. 

The  !mainf  routine  then  uses  the  pressure  functions  and  derivatives 
to  calculate  a  partial  pressure  (Equation  2.33)  for  that  wave  number  in¬ 
crement.  Each  partial  pressure  is  summed  (integrated)  so  that,  when 
the  last  wave  number  has  been  processed,  the  result  is  the  total  pres- 
sure. 

As  in  EXACT,  FFP  computes  the  propagation  loss  for  each  range. 
Implementation  of  the  FFT  would  make  it  unnecessary  to  step  through 
ranges  as  must  be  done  with  EXACT;  propagation  loss  would  be  calculated 
for  as  many  ranges  as  there  are  FFT  points. 
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IV.  TESTS  RUN  ON  THE  PROGRAM 


Because  it  was  difficult  to  gauge  the  correctness  of  the  program  or, 
more  to  the  point,  the  veracity  of  this  author's  ideas,  it  was  considered 
prudent  to  check  the  results  at  each  stage  of  development.  Although  the 
process  was  awkward  and  time-consuming  it  proved  necessary  and  resulted 
in  some  unexpected  benefits;  the  results  of  some  of  the  tests  provided 
graphical  insight  into  the  phenomenon  of  sound  propagation.  By  starting 
with  isospeed  conditions,  verifications  were  simplified.  Solutions  for 
horizontal  wave  numbers  which  resulted  from  the  computer  programs  were 
checked  against  ’correct  solutions  from  simple  formulae  for  the  non¬ 
absorption,  rigid  bottom  system,  as  well  as  for  the  system  that  included 
absorption  in  the  water.  Once  the  impedance  bottom  was  introduced,  an 
analytic  approach  was  required  to  decide  if  corrections  were  applied 
properly.  A  simple  expedient  to  check  on  the  program  at  any  point  in¬ 
volved  plotting  propagation  loss  against  horizontal  range  for  pairs  of 
modes  and  observing  the  interference  distance  between  nulls.  This  dis¬ 
tance,  when  compared  to  the  theoretical  distance  would  highlight  an 
error  if  there  were  an  anomaly.  Another  check  involved  the  observation 
that,  for  a  shallow  channel  such  as  that  used  in  the  model,  losses  are 
generally  spherical  at  short  ranges  and  cylindrical  for  far  ranges.  The 
smoothed  propagation  loss  curves  could  be  expected  to  lie  roughly  along 
curves  predicted  on  this  basis. 

The  first  tests  involved  studying  the  basic  building  blocks  of  the 
program.  As  the  data  were  read  from  input  files,  they  were  printed  into 
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a  new  file  to  ensure  that  correct  information  was  being  utilized. 
Another  simple  check  involved  printing  a  matrix  of  the  wave  numbers  as 
they  were  computed.  To  illustrate  the  importance  of  such  a  basic  check, 
it  is  worth  remarking  that  the  half-increment  method  required  that  the 
wave  number  be  calculated  for  each  half- increment  of  depth,  including  a 
depth  of  zero.  Since  the  variation  of  Fortran  used  does  not  support  a 
matrix  with  a  zeroth  element,  a  printout  proved  very  useful  in  visualiz¬ 
ing  the  situation  and  pin-pointing  a  subtle  error.  In  many  cases  it  was 
only  by  printing  out  all  of  the  variables  after  each  step  that  errors 
could  be  identified  and  remedied. 

A.  RIGID  BOTTOM  WITH  NO  ABSORPTION 

Subroutine  'eigens'  was  first  checked  for  correctness  by  comparing 
its  constant  speed  solution,  for  each  mode,  to  a  solution  which  was 
known  to  be  accurate.  This  accurate  solution  was  based  on  the  concept 
that,  for  a  constant  sound  speed,  with  rigid  bottom  and  no  absorption, 
the  pressure  function  would  be  sinusoidal,  satisfying  both  boundary  con¬ 
ditions.  This  means  that  the  horizontal  wave  number  for  each  mode  (n 
always  odd)  can  be  defined  in  terms  of  the  wave  number,  k,  and  the  water 
depth,  H: 


^  =  k2  -  (^)2.  4.1 

Because  the  subroutine  could  only  be  run  in  single  precision,  high  reso¬ 
lution  was  not  expected.  Table  4.1  shows  that  resolution  increased  lin¬ 
early  with  the  number  of  depth  increments.  For  example,  for  mode  2,  the 
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TABLE  4.1.  SUBROUTINE  ' EIGENS  '  -  COMPUTED  SQUARES  OF  WAVE  NUMBERS  AND 

ERRORS  RESULTING  FROM  VARIATION  OF  NUMBER  OF  DEPTH  INCREMENTS. 


MODE  NUMBER 


1 

2 

3 

*4 

NUM  OF  VALUE 

INCRMT  ERROR 

VALUE 

ERROR 

VALUE 

ERROR 

VALUE 

ERROR 

25 

0 .2144624294 
-0.0000906215 

0 . 1959892982 
-0.0009006298 

0.1525131932 

-0.0032877883 

0.0237859429 

-0.0614532906 

50 

0.2144174816 

-0.0000456736 

0.1955415049 

-0.0004528365 

0. 1508744489 
-0 . C016490440 

-0 . 0209414381 
-0.0167259096 

100 

0.2143947361 

-0.0000229281 

0.1953157087 

-0.0002270403 

0.1500510807 

-0.0008256757 

-0.0305268121 

-0.0071405356 

200 

0.2143832949 

-0.0000114870 

0.1952023429 

-0.0001136746 

0.1496385139 

-0.0004131089 

-0.0342955159 

-0.0033718318 

400 

0.2143775572 

-0.0000057492 

0. 1951455441 
-0.0000568758 

0 . 1494320249 
-0.0002066200 

-0.0360238122 

-0.0016435355 

800 

0.2143746840 

-0.0000028760 

0.1951171158 

-0.0000284475 

0. 14932S7310 
-0.0001033261 

-0.0368554579 
-0 . 0008118S9S 

1600 

0.2143732463 

-0.0000014384 

0. 1951028945 
-0.0000142262 

0 . 1492770720 
-0.0000516671 

-0.0372637914 

-0.0004035563 

3200 

0.2143725272 

-0.0000007193 

0. 1950957821 
-0. 0C00071137 

0. -1492512396 
-0.0000258346 

-0.0374661565 

-0.0002011912 

EXAC 

T  0.2143718079 

0.1950886633 

0.1492254049 

-0.0376673477 
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error  using  25  increments  was  ,000900,  and  the  error  was  halved  by  in¬ 
creasing  the  number  of  increments  to  50.  This  trend  continued  through 
the  run  that  used  3200  increments  and  resulted  in  an  error  of  0.000007. 
It  was  decided  that  an  error  of  0.00002,  which  resulted  from  the  use  of 
100  increments,  was  acceptable.  The  value  of  each  horizontal  wave 
number  determined  in  subroutine  'eigens'  was  intended  only  as  a  starting 
point  for  the  more  accurate  subroutines  f  modes'1  T  and  fmodes2f.  The 
accuracy,  therefore,  has  significance  only  in  that,  if  two  modes  are 

closer  together  than  single  precision  accuracy  limitations,  there  would 
be  a  chance  that  a  mode  will  be  missed  completely. 

Subroutine  Teigens!  was  also  checked  using  sound  speed  profiles  that 
varied  linearly,  both  positively  and  negatively,  with  depth.  Since  no 
clear-cut  means  of  assessing  the  error  was  available,  no  conclusions 
could  be  drawn,  and  results  are  published  for  only  the  positive  gradient 
(Table  4.2a).  In  addition,  by  using  a  sound  speed  profile  that  produced 
two  strong  ducts,  it  was  possible  to  show  that  Teigensf  was  capable  of 

distinguishing  between  two  horizontal  wave  numbers  which  are  very  close 
together  (Table  4.2b).  The  solutions  are  not  precise  but  no  modes  are 
missed  and  the  necessary  information  is  provided  to  the  subsequent  sub¬ 
routines  where  the  values  are  calculated  to  the  required  precision. 

Next,  normalized  modal  pressures  at  a  large  number  of  depth  incre¬ 
ments  were  calculated  using  tmodes2f.  These  data  were  combined  in 
Figure  4.1  and  show,  for  each  mode,  the  relative  amplitude  of  the  modal 
pressure  function.  The  profiles  for  all  three  real  modes  as  well  as  the 
first  evanescent  mode  are  plotted  together  and  care  must  be  taken  to 
note  that  they  represent  the  modal  characteris tics ,  not  the  actual 
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TABLE  4.2.  SUBROUTINE  ' E1GENS '  -  SQUARES  OF  MODAL  HORIZONTAL 

WAVE  NUMBERS  FOR  A.  POSITIVE  GRADIENT,  AND  B.  DUAL 
DUCT  FOR  50  HZ  SOURCE  AND  500  METERS  DEPTH. 


MODE 

POSITIVE 

.DOUBLE 

NUM 

GRADIENT 

DUCT 

i 

0.2023466395 

0 . 2094159634 

2 

0.2017806065 

0.2092275332 

3 

0 . 2015439631 

0.2088501530 

4 

0.2009845610 

0 . 20528282  12 

5 

0.2001803173 

0.2075239702 

6 

0.1991866971 

0 . 2055714750 

7 

0. 1980052050 

0 . 2054225531 

8 

0.1966276551 

0 . 2040741S97 

9 

0.1950384284 

0.2025220925 

10 

0. 1932226702 

0 . 2007615530 

11 

0.1911713327 

0. 1987573347 

12 

0.1888803823 

0. 1S65925S55 

13 

0.1863476142 

0 . 1941702449 

14 

0. 1835693580 

0.1915113532 

15 

0.1805375878 

0 . 1856050595 

16 

0..  1772381208 

0.1854427397 

17 

0. 1736509107 

0 . 1320079592 

13 

0.1697523042 

0 . 1782860324 

19 

0. 1655174308 

0 . 1742535547 

20 

0.1609207800 

0 . 1599023223 

21 

0 . 1559344080 

0.1551959621 

22 

0.1505243249 

0. 1501C25096 

23 

0.1446458558 

0 . 1545394076 

24 

0.1382386306 

0 . 1436057627 

25 

0.1312213613 

0. 1420936151 

26 

0. 1234851760 

0. 1349764974 

27 

0.1148818518 

0 . 1271528591 

28 

0.1051995705 

0 . 1184323554 

29 

0.0941103985 

0 . 1037542301 

30 

0.0810451735 

0. 0976845549 

31 

0.0648246885 

0.0347116344 

32 

0.0419336933 

0 . 0687321079 

33 

-0.0276373313 

0 . 0459792092 

34 

0 

-0 . 0198930514 
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pressure  contributions.  The  figure  shows  that,  for  rigid  bottom,  no  absorp¬ 
tion,  iso-speed  conditions,  the  profiles  meet  the  boundary  conditions  at 
the  top,  where  pressure  is  zero,  and  at  the  bottom,  where  the  derivative 
is  zero.  The  same  test  procedure  was  repeated  using  subroutine  'modes!' 
and  identical  results  were  obtained. 

To  illustrate  the  relative  degree  to  which  each  mode  is  stimulated, 
depending  on  the  depth  of  receiver  and  transmitter,  a  trial  run  was  made 
with  a  transmitter  at  depth  13  meters  and  receiver  at  various  depths. 
Figure  4.2  depicts  the  relative  strength  of  each  mode  at  each  depth. 

The  second  mode  is  most  strongly  stimulated  and  mode  1  is  least 
strongly  stimulated,  as  could  be  anticipated  from  Figure  4.1. 

Routine  FFP  does  not  use  exact  normal  mode  horizontal  wave  numbers, 
but  boundary  conditions  must  remain  consistent  and  sinusoidal  variations 
still  occur  at  a  rate  defined  by  the  vertical  wave  number,  kz,  which,  in 
turn,  is  a  function  of  the  horizontal  wave  number,  £.  The  results  of 
subroutine  ' modes  1'  were  compared  to  the  values  of  sin(kzz>)  and  the 
results  of  'modes2'  were  compared  to  the  values  of  cos(kz(H-z<) ,  where  H 
is  the  water  depth  and  z>  and  z<  are  the  depths  of  the  receiver  and 
transmitter.  The  comparisons  showed  that  the  subroutines  worked 
properly . 

At  this  point,  the  solution  using  calculated  horizontal  wave  numbers 
in  EXACT  was  examined  for  accuracy.  It  was  decided  that  the  accuracy 
criterion  would  be  based  upon  the  worst  case  which  allowed  for  all 
errors  to  be  cumulative.  It  can  be  seen  in  Equation  2.31  that  an  error 
in  £n  will  exhibit  itself  directly  in  the  square  root  and  in  the  exponen¬ 
tial.  A  small  error  causes  only  a  small  error  in  the  square  root  but  it 
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DEPTH  (METERS) 

50.0  46.0  40.0  36.0  30.0  26.0  20.0  15.0  10.0 


Figure  4.1.  Modal  Pressure  Function  for  Three  Real  Modes 
and  Evanscent  Mode,  Individually  Normalized. 
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DEPTH  (METERS) 

50.0  45.0  40.0  35.0  30.0  26.0  20.0  15.0  10.0 


Figure  4.2. 


Pressure  Contributions  Normalized 
of  the  Most  Strongly  Excited  Mode 


to  the  Maximum  Value 
for  Source  at  13  Meters. 
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is  critical  in  the  phase  component.  It  was  categorically  decided  that 
the  maximum  allowable  phase  error  for  each  mode  would  be  0.1  radians, 
and  that  the  maximum  range  would  be  100,000  meters.  This  limits  the 
maximum  error  in  £n  to  10“6  for  all  wave  numbers.. 

In  assessing  possible  real  errors,  a  vector  plot  was  made  for  each 
of  several  solutions.  A  constant  adjustment  was  then  introduced  to  each 
of  the  horizontal  wave  numbers,  giving  them  an  artificial  error.  It  can 
be  seen  in  Figure  4.3  that  for  short  ranges,  a  large  error  is  admissible 
and,  even  at  50,000  meters,  the  limit  imposed  by  the  accuracy  criterion 
is  well  within  the  allowable  error. 

Vector  plots  were  also  used  to  illustrate  how  a  small  change  of 
range  can  cause  a  dramatic  change  in  pressure.  Figure  4.4a  shows  the 
change  of  propagation  loss  from  830  to  840  meters.  Most  of  the  change 
is  involved  with  phase  and  the  amplitude  change  is  small.  This  can  be 
seen  best  in  Figure  4.4b  where  the  range  change  is  only  one  meter.  It 
should  be  noted  too,  that  these  errors  have  been  made  cumulative 
whereas  it  is  unlikely  that  the  errors  would  ever  be  that  way.  One 
would  therefore  expect  the  total  error  to  be  less  than  that  illus¬ 
trated. 

A  plot  of  loss  with  range  for  a  rigid  bottom,  Figure  4.5,  showed  a 
general  20  log  R  increase  in  propagation  loss  close  to  the  source  and  a 

t 

change  of  10  log  R  as  the  distance  increased.  Superimposed  on  this  is  a 
strong  interference  pattern  from  the  three  real  modes.  Sudden  changes 
of  pressure  with  changing  range,  as  were  noted  in  Figure  4.4,  are  seen  as 
rapid  fluctuations  or  propagation  loss  in  this  plot.  A  decisive  check  on 
the  general  reliablity  of  the  program  to  this  point  involved  excluding 
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Figure  4.' 


o 


Pressure  Components  for  Three  Real  Modes  Illustrating 
Cumulative  Error  from  Errors  in  Horizontal  Wave  Number 
at  50  Meters  (Upper)  and  50  Km  (Lower). 
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IMAGINARY  IMAGINARY 

0.0  1  0  2.0  3.0  4  0  -3.0  *20  -1.0  0.0 
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Components  of  Pressure  Illustrating  Sudd 
Changes  of  Pressure  with  Range. 
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Figure  4 . 4 . 
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Figure  4.5.  Propagation  Loss  for  Source  at  10  Meter  Depth 
and  Receiver  at  37  Meters  in  50  Meter  Ocean. 
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all  but  two  real  modes  and  plotting  the  resulting  propagation  loss 
curves.  For  each  pair  of  modes,  the  interference  distance  can  be  shown 
to  be  related  to  the  horizontal  wave  numbers  by: 

R  =  i,.2 

5m  5n 


where: 

R  is  the  horizontal  distance  between  nulls, 

and  £ m  and  £n  are  the  real  horizontal  wave  numbers. 

Comparison  was  made  between  the  calculated  interference  distance  and 
the  observed  distance  for  all  three  interference  combinations  and  proved 
to  be  correct  in  all  cases.  Figure  4.6  shows  the  result  when  the  first 
and  third  modes  were  used.  The  cycle  distance  was  calculated  to  be  96.5 
meters,  which  agrees  with  the  distance  from  the  plot. 


B.  RIGID  BOTTOM  WITH  ABSORPTION 

Absorption  due  to  water  (low  frequency),  although  normally  very 
small,  was  made  artificially  high  to  cheek  on  the  mechanics  of  the 
program.  With  absorption  set  to  an  arbitrary  value  of  0.0005 

nepers/meter ,  it  was  noted  that  the  imaginary  component  of  the  wave 
number  increased  with  the  increasing  grazing  angle  associated  with  higher 
modes:  0.000505  for  mode  1,  0.000555  for  mode  2,  and  0.000727  for  mode 

f 

3. 

To  ensure  that  output  losses  were  consistent  with  those  input,  other 
propagation  loss  curves  were  drawn  using  a  variety  of  rates  of  absorp¬ 
tion.  Figure  4.7  shows  that  the  difference  in  loss  at  any  given  range  is 
8.7  times  the  difference  of  absorption  rates  (as  expected  from  Equation 
2.42).  For  example,  a  change  from  0.0000001  nepers/meter  to  0.0001 
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Figure  4.6.  Interference  Pattern  Between  Modes  1  and  3 
for  Source  at  10  Meters  and  Receiver  at 
37  Meters  in  50  Meters  of  Water  for  50  Hz. 
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Figure  4.7. 


Propagation  Loss  for  Absorption  Rate  0.000001  Nepers/Meter 
(Upper)  and  0.0001  Nepers/Meter  (Lower). 
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Figure  A. 8.  Propagation  Loss  as  a  Function  of  Depth  for  Sources  at: 

a.  1  Meter,  b.  5  Meters,  c.  15  Meters,  d.  25  Meters, 
e.  40  Meters,  and  f.  45  Meters. 
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nepers/meter  causes  a  change  of  propagation  loss,  at  10  kilometers,  of 
9.0  dB.  One  would  have  expected  a  change  of  1.0  nepers  or  8.7  dB. 

A  series  of  runs  was  made  to  illustrate  how  pressure  varies  with 
receiver  depth  at  a  set  range  for  a  variety  of  source  depths.  From  the 
graphs  (Figure  4.8)  it  is  important  to  note  that  for  a  source  or  a 
receiver  at  the  surface,  the  pressure  will  be  zero  '  and  the  propagation 
loss  infinite  because  none  of  the  modes  is  stimulated.  This  does  not 
change  when  an  impedance  bottom  is  introduced,  but  it  would  change  if  a 
real  surface  were  considered. 

Using  these  same  depth  profiles,  it  was  possible  to  further  verify 
the  program  by  checking  for  reciprocity.  The  propagation  loss  is,  for 
example,  the  same  for  a  combination  of  a  source  at  5  meters  and  receiver 
at ’25  meters  as  for  a  source  at  25  meters  and  receiver  at  5  meters. 

The  final  set  of  tests  for  a  rigid  bottom  involved  the  use  of  the 
FFP.  At  the  point  of  testing,  the  Fourier  transform  ( FFT)  had  not  yet 
been  introduced  so  pressure  factors  were  being  calculated  for  individual 
ranges.  Figure  4.9  shows  that,  for  a  zero  absorption  loss,  the  only  hor¬ 
izontal  wave  numbers  that  contribute  to  the  total  are  at  the  exact  mode 
values.  Small  changes  of  the  horizontal  wave  number  cause  serious 
changes  in  the  value  close  to  the  exact  mode  value. 

Introduction  of  absorption  causes  the  energy  to  be  dispersed  across 
the  spectrum.  There  are  small  contributions  from  across  the  wave 

number  spectrum  and  flattening  at  the  modes.  Higher  absorption  causes 

greater  dispersion  of  the  pressure*.  The  integral,  with  absorption, 
can  be  seen  to  fluctuate  in  Figure  4.10,  staying  close  to  zero  until  the 

first  mode  is  approached.  The  integral  increases  quickly,  but  not 
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PRESSURE 


Figure  4.9. 


Amplitude  of  Sound  Pressure  Factor  for  Individual 
Wave  Numbers  with  No  Absorption. 
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PRESSURE  INTEGRAL 

0.0  6.0  10.0  15.0  20.0  25.0  30.0  36.0 


Figure  4.10.  Integrated  Pressure  with  Absorption  Rate 

0.00005  Nepers/Meter  at  Range  10000  Meters, 
Tx/Rx  at  10/37  Meters,  50  Meter  Ocean,  50  Hz. 
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abruptly  when  a  mode  is  reached.  A  sharp  decline  is  notable  when  the 
1 80  degree  phase  change  occurs  at  the  mode.  The  final  value  represents 
the  final  pressure  at  the  receiver. 

Various  checks  were  made  to  ensure  proper  functioning  of  the  FFP 
program  for  a  given  range.  By  changing  the  number  of  wave  number  incre¬ 
ments,  it  was  found  that  the  maximum  amplitude  for  each  mode  changed 
(Figure  4.11)  but  that  the  integrated  pressure  was  identical  when  the 
number  of  increments  was  changed  from  1024  to  2048.  Pressures  agreed 
with  one  another  to  the  seventh  decimal  and  fifth  significant  digit. 
Given  more  time,  it  would  have  been  interesting  to  see  how  few  samples 
could  be  taken,  for  a  single  range,  before  significant  errors  were  intro¬ 
duced. 

Figure  4.12  shows  the  effect  of  changing  absorption  rates.  However, 
the  change  of  pressure  at  a  given  range  does  not  correlate  with  the 
change  of  absorption  rate.  This  indicated  that  the  FFP  routine  was  in¬ 
correct  in  some  way  and  further  checks  were  required. 

The  integrated  pressures  from  this  program  were  then  compared  to 
the  values  obtained  from  the  program  EXACT  for  different  absorption 
rates  and  ranges.  As  can  be  seen  in  Table  4.3,  agreement  was  poor. 
Because  of  time  constraints,  it  was  decided  to  abandon  the  FFP  and  con¬ 
centrate  on  solutions  using  the  EXACT  method. 

NOTE:  A  great  deal  of  time  had  been  spent  on  FFP  trying  to  find  a  solu¬ 
tion  for  a  fixed  range,  a  prerequisite  to  implementation  of  the  'FFP' 
subroutine.  At  the  time  that  work  with  the  FFT  was  first  suspended,  it 
gave  inconsistent  results  that  did  not  meet  any  of  the  testing  criteria. 
A  small  change  of  increment  size,  from  2047  to  2048  increments,  caused 
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TABLE  4.3.  CALCULATED  PRESSURE  FOR  PROGRAM  FFP  AND  EXACT. 


RANGE 

ABSORPTION 

PRESSURE 

FACTOR 

(METERS) 

(NEPERS/METER) 

FFP 

EXACT 

1  000 

0.0000001 

0.001  72867 

0.00139872 

1  000 

0.00001 

0.001  281  78 

0.001 38953 

2000 

0.0000001 

0.00001 486 

0.00139504 

2000 

0.00001 

0.0001 0503 

0.001  381  01 
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PRESSURE  FACTOR 


HORIZONTAL  WAVE  NUMBER 


Figure  4.11.  Amplitude  of  Sound  Pressure  for  Absorption 
0.0005  Nepers/Meter  Using:  1024  Wavenumber 
Samples  (Dashed  Curve)  and  2048  Samples  (Solid.) 
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PRESSURE  INTEGRAL 

0.0  6.0  10.0  15.0  20.0  25.0  30.0  35.0 


°  =  .00001 
o  =  .0000001 
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HORIZONTAL  WAVE  NUMBER 


Figure  4.12.  Integrated  Pressure  for  Absorption  Rate 

a.  0.00001  Nepers/Meter  (Lower)  and 

b.  0.0000001  Nepers/Meter  (Upper)  Range  1000  Meters. 
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outlandish  changes  in  the  integrated  pressure.  The  changes  of  pressure 
resulting  from  a  changed  absorption  rate  did  not  correspond,  even  rem¬ 
otely,  with  the  changes  expected.  At  that  point,  priorities  were  shifted 
and  it  was  decided  to  concentrate  on  solving  the  'EXACT'  program  prob¬ 
lems  at  the  expense  of  the  'FFP'.  At  a  later  date,  too  late  to  be  prac¬ 
tically  helpful,  some  of  the  FFP  problems  were  rectified  and  the  above 
mentioned  inconsistencies  were  remedied.  It-  was  possible  to  relate 
propagation  loss  to  absorption  rate  and  it  was  shown  that  the  number  of 
samples  has  a  negligible  effect  on  the  final  pressure.  However,  the 
discrepancy  between  the  EXACT  solution  and  the  FFP  solutions  still  exist 
and  the  FFP  is  incomplete. 
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C.  IMPEDANCE  BOTTOM 


The  introduction  of  an  impedance  bottom  to  replace  the  perfectly 
reflecting  rigid  bottom  made  it  possible  to  better  predict  the  propaga¬ 
tion  loss  experienced  in  reality.  This  was  done  by  defining  the  bottom 
in  terms  of  its  density  and  the  sound  speed  in  the  bottom  at  the  inter¬ 
face.  After  continued  failure  of  the  program  to  solve  for  modes  whose 
grazing  angles  were  far  below  the  critical  angle  and  the  eventual  reso¬ 
lution  of  that  problem,  a  series  of  tests  was  devised  to  try  to  optimize 
the  program.  These  tests  also  provided  a  convenient  means  of  illustrat¬ 
ing  how  various  angles  of  incidence  associated  with  the  different  modes 
resulted  in  the  varying  phase  changes  and  amplitudes  of  reflected 
energy. 

As  with  the  rigid  bottom  case,  it  was  possible  to  plot  the  modal 
pressure  as  a  function  of  depth,  from  the  reflecting  surface  to  the  imp¬ 
edance  bottom.  Figure  4.13  shows  how  the  pressure  reaches  a  maximum 
well  above  the  bottom  for  mode  1.  The  grazing  angle,  in  this  case,  is 
less  than  the  critical  angle  and  illustrates  the  phase  advance  at  the 
bottom.  The  test  was  repeated  for  strong  positive  and  negative  gradi¬ 
ents  but  showed  little  change. 

The  first  set  of  tests  to  check  on  final  accuracy  involved  using  a 
critical  angle  of  0.3  radians  and  was  intended  to  find  the  optimal 
number  of  increments.  It  had  been  determined  earlier  that  accuracy  to 
the  ninth  decimal,  or  better,  was  obtainable  (iso-speed)  for  all  modes 
using  200  depth  incements/mode.  At  100,  the  accuracy  dropped  to  the 
seventh  decimal  and  at  50,  answers  were  accurate  to  only  the  fourth 
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DEPTH  (METERS) 

60.0  46. 0  40.0  36.0  30.0  36.0  20.0  16.0  10.0  6.0  0  0 


} - 1 -  |  y  w  'V 

0.00  C.25  0.50  0.75  l.CC 

RELATIVE  PRESSURE  FUNCTION 


Figure  4.13.  Relative  Pressure  Function  for  3  Real  Modes  for  Critical 

Angle  0.3  Radians,  Isospeed  Water  Conditions,  30  Hz  Source. 
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decimal.  Table  4.4  shows  horizontal  wave  numbers  that  were  calculated 


using  100,  200,  400  and  600  depth  increments/mode  in  combination  with 
correction  factors  of  1,  2,  5  and  10.  The  correction  factor  was  divided 
into  the  mode  wave  number  correction  term  of  subroutines  'modesl'  and 
'modes2'to  limit  the  size  correction  and  thus  prevent  instability  . 
These  correction  factors  were  found  necessary  for  modes  whose  grazing 
angles  were  less  than  the  critical  angle  (in  this  case  only  mode  1). 
Without  a  correction  factor,  the  program  would  fail  to  converge  to  an 
acceptable  solution  but  would  become  divergent. 

A  major  problem  in  deciding  the  best  choice  of  variables  came  about 
because  there  was  no  correct  answer  against  which  to  appraise  the 
values.  It  was  reasoned  that  800  increments  per  mode  and  a  correction 
factor  of  10  would  give  the  most  accurate  answer  and  the  problem  nar¬ 
rowed  down  to  finding  a  combination  that  gave  an  acceptable  accuracy 

with  a  reasonable  number  of  iterations.  Having  decided  that  200  incre¬ 
ments/mode  was  sufficient  and  necessary  for  the  rigid  bottom  case,  it 
became  a  question  of  how  large  the  correction  factor  could  be.  The 

choice  was  further  narrowed  to  a  correction  factor  of  2  with  an  error  of 
9.8  x  10-8  and  40  iterations  required  and  a  factor  of  5  with  an  error  of 
3  x  10~8  and  89  iterations.  Because  the  error  was  borderline  for  the 

correction  factor  2,  the  final  choice  was  200  intervals/ mode  with  a  cor¬ 
rection  factor  of  5. 

Using  the  inputs  as  decided  above,  an  effort  was  made  to  correlate 
actual  accuracy  with  that  specified  as  a  requirement.  Results  for 

trials  on  all  modes  using  accuracy  criteria  from  10~5  to  10~13  follow  in 
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TABLE 


.4.  PROGRAM  'EXACT'  -  REAL  AND  IMAGINARY  COMPONENTS  OF  £2 
FOR  MODES  1  TO  4  RESULTING  FROM  VARIATION  OF  NUMBER 
OF  DEPTH  INCREMENTS  AND  'CORRECTION'  FACTORS. 


COR  i:;c  ITER 


1 

100 

90 

2 

100 

38 

5 

100 

39 

10 

100 

185 

1 

200 

74 

2 

200 

40 

5 

200 

£9 

10 

200 

195 

x 

400 

67 

2 

400 

40 

5 

400 

89 

10 

o 

o 

-i* 

186 

1 

800 

72 

2 

800 

39 

5 

800 

90 

10 

800 

187 

EXACT 

(RIGID  BOTT) 


0.2I2512S34 

0.000450164 

0.212512965 

0.000450172 

0.212513034 

0.000450172 

0.212513050 

0.000450170 

0.212512355 

0.000450155 

0 . 212512368 
0.0C0450171 

0.212513035 
0 . 00C450 172 

0.212512063 
0 . 0C0450 170 

0.212512555 
0 . 0C0450155 

0.212512959 

0.0004501"! 

0.212513023 
0 . 000450 1"2 

C. 212512065 
0.000450170 

0.212512S85 

0.000450165 

0.212512959 

0.000450171 

0.212513039 

0.000450171 

0.212513066 

0.000450170 

0.214371320 

0.000505340 


0 . 195507949 
-0 . C0333983 1 

0. 1955C7949 
-0 . C0333933 1 

0. 195507949 
—  0 . 00333983 1 

0. 195507949 
-0.003339331 

0. 195507947 
-0 . CC3339330 

0. 1955C7947 
-0.003339330 

C.  195507947 
-0.003239830 

0. 195507947 
-0.003339830 

0. 195507946 
-0 . 003339330 

0. 195507946 
-0.003339830 

0.  195507946 
-C . 003339330 

0 . 195507946 
-0.003339820 

0.  195507946 
-0 . 003339830 

0. 195507945 
-0 . 0033  2983  0 

0. 195507946 
-0.003339830 

0.195507946 

-0.003339830 

0.195088818 

0.000555269 


C.  150035436 
-0.010172100 

0 . 15003  543  6 
-0 . 01C172100 

0. 150035436 
-0:010172100 

0. 15C03  5436 
-0 . C1C1721C0 

0.  15C03  5443 
-0.0 1C  172094 

0. 150035443 
-C. 010172094 

0 . 150C35443 
-0.010172094 

0. 150035443 
-0.010172094 

0. 150035443 
-0 . 0101 72093 

0. 150035443 
-0.010172093 

0. 150035443 
-0.010172093 

0. 150035443 
-0.010172093 

0. 150035443 
-0.010172092 

0. 150035443 
-0. 010172092 

0. 150035443 
-0.010172092 

0. 150035443 
-0.010172092 

0. 149226333 
0.000725949 


0.002967414 

0.037779947 

0.002367414 

0.037779947 

0. C02SS7414 
0.037776947 

0.002867414 
C. 037779947 

0 . 002867435 
0 . 02  77  79660 

0 . 002367435 
0 . 0377" 9 660 

0  . CG255"43 5 
0. 03"""9660 

0 . CC2367435 
C. 037779650 

0.002867437 
C. 037779642 

C . C02867437 
0 . 037779642 

C . 002  36743  7 
C. 037779642 

C  002357*- 3 n 
0 . 03 ""79642 

0.002867437 
C . C3  7779541 

C . 00236743  7 
0.037779641 

0. 00286~437 
0.037779541 

0.002367437 
0 . 037779541 

0.002867437 
0 . 03  77"9540 
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Table  4.5.  Assuming  that  the  answers  obtained  using  10~13  had  the  least 
error  and  could  then  serve  as  a  standard,  errors  were  calculated  as  the 
difference  between  the  standard  and  each  answer.  It  was  found,  in  all 
cases,  that  the  actual  errors  were  slightly  less  than  those  specified  by 
the  accuracy  criteria.  This  was  to  be  expected  since  corrections  were 
made  to  make  the  error  less  than  the  designated  amount. 

The  final  check  involved  the  problems  associated  with  low  modes 
whose  grazing  angles  were  less  than  the  critical  angle.  A  series  of 
calculations  was  made  to  find  the  number  of  iterations  required  for 
different  critical  angles  (different  sound  speed  ratios)  using  a  100 
meter  deep,  isovelocity  ocean  that  supported  seven  real  modes.  It  was 
also  intended  to  uncover  further  problems  associated  with  large  critical 
angles.  At  the  same  time,  the  squared  values  of  horizontal  wave  numbers 
for  the  seven  modes  were  calculated  for  each  critical  angle  and  are 
shown  in  Table  4.6.  For  angles  close  to  the  grazing  angle  for  the  first 
mode,  0.1  radians,  the  program  failed  using  both  subroutines  'modesl'  and 
'modes2'  during  these  runs.  This  was  unusual  because  previous  data  had 
been  collected  for  almost  identical  conditions. 

The  fraction  of  energy  reflected  and  the  phase  change  represented  by 
the  Rayleigh  coefficient,  were  calculated  for  each  mode  for  each  criti¬ 
cal  angle  and  are  shown  in  Figure  4.14  and  4.15  (upper).  The  graphs  show 
that  the  loss  per  bounce  is  zero  when  the  grazing  angle  is  less  than  the 
critical  angle.  Mode  1  has  the  lowest  grazing  angle  and  it  can  be  seen 
that  the  critical  angle  must  be  very  low  if  the  low  modes  are  to  suffer 
any  attenuation  at  all.  Higher  modes  suffer  higher  losses.  To  show  the 
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TABLE  4.5.  PROGRAMS  ' EXACT'  -  REAL  AND  IMAGINARY  COMPONENTS  OF  ^ 
AND  NUMBER  OF  ITERATIONS  TO  MEET  REQUISITE  ACCURACY. 


MODE  NUMBER 


1 

2 

3 

ACC 

REAL 

I  MAG 

ITERATIONS 

REAL 

I  MAG 

ITERATIONS 

RE.AL 

I  MAG 

ITERATIONS 

REAL 

I  MAG 

ITERATIONS 

E-5 

0. 21251957S303 
0.000446662283 
18 

0 . 195507093057 
0.0033392S5283 

6 

0. 150034762633 

0 .0101639743 14 

5 

0 . 0028674’71359 
0.037779636753 
3 

E-7 

0.212512933243 

0.000450155917 

0. 195507946339 
0.003339829696 
£ 

0. 150035442823 
0.010172092319 

7 

0.002867437052 

0.037779642325 

4 

E-a 

0. 212512S91707 
0.000450175339 
38 

0.195507923776 

0.003339803692 

9 

0. 150025443413 
0.010172073562 

8 

0.0C2867437052 

0.037779642325 

4 

E-9 

0.212512888681 

0.000450176659 

44 

0.195507928795 

0.003339804619 

11 

0. 150035441633 
0.010172073282 

9 

0.002S67437109 

0,037779642316 

5 

E  - 10 

0.212512888032 

0.000450176927 

51 

0. 195507928955 

0. 003339304648 

12 

0. 150035441537 
0.010172073444 
10 

0.002367 43 7109 
0.037779642316 
5 

E-  12 

0.212512387979 

0.000450176947 

58 

0. 195507928962 
0.003339804619 
14 

0. 150035441601 
0.010172073451 

11 

0.002367437109 

0.037779642316 

5 

E-  13 

0.212512887974 

0.000450176949 

70 

0. 195507928961 
0.003339804620 

16 

0. 150035441602 
0.010172073449 

13 

0.002867437109 

0.037779642316 

6 
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PHASE  CHANGE  (RADIANS) 


CRITICAL  ANGLE  (RADIANS) 


Figure  A. 14.  Phase  Change/Bounce  for  3  Real  Modes  as  a  Function 
of  Critical  Angle. 
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TABLE  4.6. 


PROGRAM  'EXACT'  -  REAL 
FOR  IMPEDANCE  BOTTOM  AS 
(0  TO  .25  RADIANS) . 


AND  IMAGINARY  COMPONENTS  OF  r,2 
A  FUNCTION  OF  CRITICAL  ANGLE ’ 


MODE  NUMBER 


CRIT 

ANGLE 

1 

2 

3 

4 

5 

6 

7 

RIG 

0.214751 

0.210105 

0. 200490 

0.165134 

0.162416 

0. 128489 

0 . 068307 

BOT 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

O.OCOOCC 

.25 

0.209522 

-0.209622 

0.20C611 

0. 185266 

0. 162595 

0 . 128830 

0.070372 

-0.000004 

-0.000004 

-0.001557 

-0.002882 

-0.004518 

-0.007198 

-C. 015854 

.20 

0.210226 

0.210226 

0 . 200600 

0.185258 

0.162589 

0. 128825 

0.070365 

-0.000471 

-0.000471 

-0.001785 

-0.003022 

-0.004612 

-0 . 007258 

-0.015682 

.  15 

0.214297 

0.210212 

0 . 200581 

0.185252 

0. 162584 

0. 128821 

0.070360 

-0.000004 

-0.000863 

-0.001949 

-0.003126 

-0.004682 

-0.007303 

-0.015902 

.  10 

UNAVAILABLE 

UNAVAILABLE 

.05 

0.213628 

0.213628 

0. 198860 

0.185216 

0.162580 

0.128792 

0.07CC67 

-0.000499 

-0.000499 

-0.000504 

-0.000591 

-0.002909 

-0.005782 

-0.014147 

o 

o 

0 . 214649 

0.210107 

0.200526 

0.185199 

0.162529 

0. 128739 

0 . 069593 

-0.000123 

-0.000712 

-0.001612 

-0.002578 

-0.004114 

-0.006535 

-0.C14427 
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extreme  modal  dependence  of  attenuation,  Figure  4.15  (lower)  displays 
the  loss  per  meter  for  each  mode  for  a  variety  of  cirtical  angles. 

These  curves  illustrate  that  the  increased  grazing  angles  associated 
with  the  higher  modes  result,  not  only  in  higher  loss  per  bounce,  but 

also  in  an  increased  number  of  bounces.  For  a  small  critical  angle,  say 
0.1  radians,  the  loss  per  bounce  is  approximately  three  times  as  great 
for  mode  2  as  for  mode  1,  but  the  loss  per  meter  is  about  ten  times  as 
great  (Figure  4.15). 

It  was  found  that  a  change  of  critical  angle  (change  of  the  sound 
speed  ratio)  changes  the  loss  greatly.  Figure  4.16  demonstrates  the 
change  of  pressure  due  to  mode  1  when  the  critical  angle  is  changed  from 
above  to  below  the  grazing  angle.  Changing  the  impedance  ratio  by  ten 
per  cent  caused  a  corresponding  increase  in  the  loss  but  did  not  change 
the  point  at  which  the  critical  angle  became  effective. 

To  study  the  bottom  loss  effects  for  grazing  angles  greater  than 
critical,  different  runs  were  made  for  critical  angles  that  allowed  a 

different  number  of  modes  to  propagate  in  the  low  loss  manner  (grazing 

angle  less  than  critical  angle).  Figures  4.17  and  4.18  show  how  much 
each  mode  contributes  toward  the  total  pressure,  depending  primarily 
upon  whether  its  grazing  angle  is  above  or  below  the  critical  angle. 
They  show  that  for  each  mode,  when  the  grazing  angle  is  less  than  the 
critical  angle,  the  propagation  loss  for  that  mode  is  identical  to  that 
for  the  rigid  bottom  case. 

As  a  last  demonstration  of  the  effects  of  the  impedance  bottom  on 
propagation,  data  were  compared  for  the  isovelocity  case  and  the  strong 
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o 


Figure  4.13*  Loss  Bounce  (Upper)  and  Loss/Meter  (Lower) 
Due  to  Bottom  Interaction. 
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DEPTH  (METERS) 

60. 0  46.0  40  0  38.0  30.0  26.0  20.0  16.0  10.0  6.0  0.0 


Figure  4.16.  Pressure  Function  for  Mode  1  With  Critical  Angle  Above 
Grazing  (Crit  =  0.30)  and  Below  Grazing  (crit  =  0.10). 
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PROPAGATION  LOSS  (Dll)  PROPAGATION  LOSS  (DB) 

120.0  110. 0  100.0  00  0  00.0  70.0  00.0  J20.0  110.0  100.0  eo.o  60.0  70.0  00.0 


DISTANCE  (METERS) 

Figure  4,17.  Propagation  Loss  for  Each  Mode  When  Critical  Angie 
is  Less  Than  Grazing  Angle  for  All  Modes  (Upper) 
and  For  Mode  1  Only  (Lower). 
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PROPAGATION  LOSS (DO)  PROPAGATION  LOSS (DO) 

100.0  90.0  ao.o  70.0  60.0  120.0  110.0  100.0  00.0  60. 0  70.0  60.0 


DISTANCE  (METERS) 


Figure  A. 18.  Propagation  Loss  Cor  Each  Mode  When  Critical  Angle 

Less  Than  Grazing  Angle  tor  Modes  1  and  2  Only  (Upper) 
and  for  Modes  1,  2,  and  3  (Lower).  . 
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positive  gradient  case  when  the  critical  angle  was  0.1  radians.  It  can 
be  seen  in  Figure  4.19  that  the  upward  refraction  casued  by  the  gradient 
lessens  loss  due  to  bottom  interactions. 
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PROPAGATION  LOSS (DB) 

100.0  90.0  00.0  70.0  00.0  GOO  40.0 


Figure  4.19.  Propagation  Loss  with  Impedance  Bottom  for 
Isovelocity  Water  Conditions  (Dashed  Lines) 
and  Strong  Positive  Gradient  (Solid  Line). 
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V.  CONCLUSION 


The  program  EXACT,  which  computes  propaation  loss  by  first  solving 
for  individual  mode  wave  numbers,  was  found  to  give  realistic  values. 
The  Fast  Field  Program,  which  was  meant  to  solve  for  loss  using  an  FFT 
was  not  completed  because  of  complications  that  proved  unsolvable  in  the 
limited  time  available. 

One  of  the  problems  with  the  solution  of  the  program  EXACT,  which 
would  also  be  a  problem  with  FFP,  is  the  restriction  placed  on  the  impe¬ 
dance  bottom.  All  of  the  sound  transmitted  into  the  bottom  is  consi¬ 
dered  lost  to  the  bottom.  Mo  allowance  is  made  for  this  energy  to  be 
refracted  upward  and  re-transmitted  into  the  water.  This  might  not  be 
entirely  unrealistic  in  some  cases,  because  bottom  absorption  is  often 
very  high.  However,  this  simplification  limits  the  practicality  of  this 
program  and  is  considered  by  the  author  to  be  its  main  weakness. 

Another  point  of  consideration  is  the  number  of  modes  involved.  The 
tests  run  usually  involved  only  three  real  modes,  two  of  which  are  lost 
to  the  bottom  after  a  short  range  when  the  critical  angle  is  0.1  radians 
(sound  speed  ratio  is  1.005:1).  An  increase  of  depth  or  frequency  will 
mean  an  increased  number  of  modes.  There  might  be  it,  400,  or  40,000 
modes  present  but  a  large  portion  of  them  may  suffer  a  large  bottom 
loss.  This  also  explains  why  trapped  modes  are  so  significant  and  why 
it  is  often  sufficient  to  consider  only  a  few  modes  in  most  analyses. 


78 


One  of  the  intended  purposes  of  the  thesis  was  to  compare  the  com¬ 
puter  processing  time  for  each  program  to  solve  for  conditions  where  a 
large  number  of  modes  were  present.  It  was  reasoned  that  an  increase 
of  frequency  and/or  depth  would  require  a  proportionate  increase  in  the 
number  of  depth  samples  used  by  the  subroutines  that  calculated  pres¬ 
sure  as  a  function  of  depth  for  each  horizontal  wave  number.  This  would 
mean  that  if  the  frequency  were  increased  from  50  to  100  Hertz  and  the 
depth  from  50  to  500  meters,  there  would  be  a  20  fold  increase  in  the 
processor  time  used  to  do  these  calculations  in  each  of  EXACT  and  FFP. 
However,  there  would  be  a  further  20  fold  increase  in  time  for  EXACT 
because  it  would  be  required  to  calculate  for  20  times  as  many  modes. 
The  FFP  would  not  require  any  similar  increase  and  would  therefore  have 
a  decided  advantage  when  a  large  number  of  modes  were  involved.  It  is 
possible  that  for  cases  with  a  small  number  of  modes,  EXACT  would  be 
faster,  but  the  information  was  not  available  for  comparison. 

A  simplification  in  EXACT  involves  the  accuracy  criterion  of  the  hor¬ 
izontal  wave  number.  It  was  decided  that  errors  due  to  accuracy  limita¬ 
tions  would  be  acceptable  if  the  wave  number  were  resolved  to  less  than 
10-6.  However,  the  subroutines  'eigens',  'modesl',  and  'modes2'  return 
the  square  of  the  wave  number.  For  the  calculated  £n  to  remain  within 
the  limit, 


^  +  |  Error  |  <  Un  +  10"6)2 
SO 

| Error  |  <  2  x  10-s  £n. 
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For  the  conditions  considered,  the  correct  accuracy  criteria  should  have 
been  marginally  more  restrictive  for  all  three  modes  considered.  Had 
any  very  small  wave  numbers  resulted,  the  allowable  accuracy  criteria 
should  have  been  much  more  restrictive.  A  more  complete  program  would 
calculate  each  wave  number  to  its  own  allowable  accuracy  limit  using 
formula  5-1 

While  it  is  recognized  that  program  results  would  in  no  way  alter 
the  facts  of  nature,  the  correspondence  between  observed  phenomena  and 
predicted  results  did  add  to  the  credibility  of  the  theory  and  program 
outlined  in  this  thesis.  Many  of  the  results,  as  illustrated  in  the 
graphs,  would  prove  useful  as  teaching  aids  to  elucidate  the  principles 
of  sound  propagation  in  the  ocean.  Without  reference  to  the  development 
of  normal  mode  theory,  it  is  possible  to  illustrate  the  consequences  of 
changing  depth,  range,  frequency,  water  depth,  and  sound  speed  profiles. 
For  example,  the  program  could  be  used  to  show -the  effects  of  surface 
duct  on  long  range  propagation. 

Finally,  it  is  realized  that  the  programming  methods  used  are  far 
less  efficient  than  many  readily  available  programs.  The  intention  of 
this  thesis  was  to  develop  an  individual  Normal  Mode  program  from  a  the¬ 
oretical  basis  without  direct  reference  to  other,  more  sophisticated  pro¬ 
grams.  A  more  complete  study  would  compare  program  results  to  actual 
observed  data  as  well  as  to  results  of  other,  more  authoritative 
predictions . 
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APPENDIX  A 


MATRIX  METHOD  FOR  FINDING  MODAL  HORIZONTAL  WAVE  NUMBERS 

The  matrix  method  (Ref. 8)  employed  in  the  program  EXACT  yields  a 
complete  set  of  eigenvalues  in  a  single  pass  utilizing  IMSL  routine 
EQRTIS.  It  is  convenient,  fast  and  has  the  advantage  that  all  modes  are 
assuredly  determined,  no  matter  how  close  together  they  are.  Accuracy 
is  not  attempted;  that  is  assured  with  the  Runge-Kutta  routines  that  are 
used  subsequently. 

The  Matrix  solution  is  a  means  of  solving  for  the  eigenvalues  of  the 
one  dimensional  differential  equation: 


where: 

4>  is  the  eigenfunction, 


kz  is  the  vertical  wave  number  / (k £*)  , 

5  is  the  horizontal  wave  number,  and 
z  is  the  depth. 

It  solves  the  equation  by  finding  the  eigenvalues  of  the  matrix  M  in 
the  following  Matrix  equation: 


mijj  =  -(Az)2  Et|i 


A. 2 


81 


This  equation  results  from  writing  the  differential 
difference  equation  at  each  of  N  depth  Az  apart.  >l> 

N  elements ,• each  element  defined  by: 


equation  as  a  finite 
a  column  matrix  with 


i*;n  =  >KnAz)  n  =  1  ,  2 - N 


A. 3 


and  M  is  a  square  symmetric  tridiagonal  matrix  which  allows  a  computation¬ 
ally  efficient  means  of  finding  the  eigenvalues. 


D,  -1  0 

-1  D2  -1 

0  -1  D3 


0 

0 


0 

0 


where: 


A. 4 


Dtr1  -1 
0  - 1  Djj 


Dn  =  2  -  (Az)2k2(nAz) 


n  =  1,  2, - ,N 


A. 5 
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APPENDIX  B 


COMPUTER  PROGRAM  EXACT 


.rtkam  .m 


C-  MAIN  PROGRAM  FUADS  UP  TO  30  PAIRS  OF  SOUND  VELOCITY  INFORMATION’ 

C*  IT  ASSIGNS  CALCULATED  VALUES  OF  WAVENUMBER  FOR  EACH  INCREMENTAL 

C*  DEPTH  AMD  CALCULATES  THE  NUMBER  OF  REAL  EIGENVALUES  AS  WELL  AS 
C*  MAX  WAVENUMBER  AND  THE  DEPTH  INCREMENT 

C*  ABSORPTION  IS  ASSUMED  -TO  EE  CONSTANT  AT  ALL  DEPTHS 

Cw  NON-COMPLEX  EIGENVALUES  ARE  CALCULATED  BY  A  MATRIX  METHOD  IN 

C*  SUBROUTINE  El GENS  AND  THESE  ARE  USED  AS  A  FIRST  ESTIMATE  IN  ONE 
C*  OF  MODES  SUBROUTINES 

C-  PRESSURES  AND  THE  WRONSKI AN  FROM  MODES  ARE  USED  TO  CALCULATE 
C*  COMPLEX  INDIVIDUAL  MODAL  PRESSURES  WHICH  ARE  SUMMED  FOR 
C*  PARTICULAR  RANGES 

C-  BOTTOM  ABSORPTION  NOT  ACCOUNTED  FOR 

C* 

COMPLEX  K(  10001) , P ( 10001 ) , CALC ( 50 ) , ARGUE  , VAL  , CORR , I NT, ?D I FF , PTM , 
1PRX , DENOM, PPRESS ( 50 ) , EIG(50) , PRESS (50) , PRES SR, IMAG , ALFA , AIG ( 50 ) , 
1 FACT, REF, COMP 

REAL  D ( 30 ) , KAY ( 30 ) , RATE ( 30 ) , C { 30 ) , E IG2 ( 50 ) , MAG ( 50 ) , KHAX, Ki , 

1KAE( 10001 ) , KAA( 10001 ) , PR(5C) , PIN (50) 

LOGICAL  TO? 

C 

TOP  =  .FALSE. 

NUM=1 

PI  =  3 . 141592654 

CKAX=0 

CMIN=200 

FREQ  =  50 

DEPTHR  =  10 

DEPTHT  =37 

RMIM  =20 

INC  =  400 

JUMP  =  INC/100 

ERRMAX  =  .00000001 

CORN  =  2 

CRITA  =  .15 

RATIO  =2.0 

12  READ  ( 3 , *  )  D(NUM) , C(NUM) 

C 

C-  EACH  SOUND  VELOCITY/DEPTH  PAIR  HAS  A  WAVENUMBER  CALCULATED 
C 

IF  (D(NUM)  . LT.  0)  GOTO  14 
KAY(NUM)  =2  *  PI  *  FREQ  /  C(NUM) 

IF  (C(NUM) .GT.CMAX)  THEN 
CMAX=C ( MUM ) 

END  IF 

IF  ( C ( MUM ) . LT . CM I N )  THEN 
KMAX=KAY ( MUM ) 

CM IN  =  C(NUM) 

END  I  F 

IF  (MUM  . EQ.  1)  GOTO  13 
C 

RATE  (MUM-1)  =  (KAY (MUM)  -  KAY (MUM-1))  /  ( D ( MUM ) -D ( MUM- 1 ) ) 
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P  i  :  .ia  ACT  FOR"' RAN 


PRINT  :03,D(NTJi-:-l ) , C ( NUN- I ) , KAY ( NUN- 1 ) , RATE ( NUN- 1 ) 

13  NUM  =  NUM  +  1 
GOTO  12 

14  BOTTOMED ( NUN- 1 ) 

N  =  I  FIX ( ( 4*BOTTOMx  FREO/CMAX  -  l)/2.0) 

WRITE  (  4 , 63  89  )  M 

6339  FORMAT  (//'  THERE  ARE  ',16,'  REAL  MODES'/) 

M  =  1  +  N 

DELZ  =  BOTTOM/ ( I NC*N ) 

PRINT  103 , D ( NUM- 1 ) ,C(NUM-1) 

WR I TE  (  4 ,  1 03 ) D ( MUM- 1 ) # C ( NUM- 1 ) , KAY ( NUM- 1 ) 

IF  (C(NUM-l) .LT.C(l) )  TOP  =  .TRUE. 

IF  (N.GT.50)  THEN 
WRITE (4, 100)  N 

100  FORMAT ( 1  TOO  MANY  MODES  FOR  THIS  PROGRAM!!  M ( REAL )  =  ’,14) 

GOTO  11 
END  IF 

x  rutilrjfXTximrifitxitjniMiiitrrrTixj  r-xrxx-xxxxxxx 

* 

*  THE  NEXT  SECTION  CALCULATES  THE  EXACT  WAVENUMBERS  FOR  THE  SU3- 

*  ROUTINES  E I GENS ,  MODES1  AND  MODES 2 

*  THE  REAL  COMPONENTS  FOR  MODES 1  AND  MODES 2  ARE  CALCULATED  FOR 

*  INC*N  +  1  DEPTHS  STARTING  WITH  DEPTH  1  AT  THE  SURFACE  AND  INC* 

-  AT  THE  BOTTOM 

k  AT  REGULAR  DEPTH  INCREMENTS  A  REAL  VALUE  OF  K  IS  TAKEN  FOR 

*  USE  IN  SUBROUTINE  E I GENS  THE  NUMBER  OF  SAMPLES  IS  ICO 


KOUNT  =  C 
NUM  =  1 

DO  18  IG  =  1,  I NC  *  M  + 1 
DP  =  (IG-1)  -  DELZ 

K  A  A  (  I G  )  =  KAY  ( NUM )  +  RATE  ( I  iUM )  *  (  DP  -  D  ( i  .‘UN  )  ) 

IT  IS  POSSIBLE  TO  INSERT  A  LOOP  TO  CHOOSE  ONLY  A  FRACTION  OF 
THE  DEPTH  SAMPLES  OF  K(  )  A  TEST  RUN  WAS  MADE  TC  COMPARE  THE 
VALUES  OBTAINED  FOR  VARIOUS  'JUMPS' 

rxx*xx:Vx*xxxxxxx~x*x:TXXxxxxxx*xxxx7r^xx*x:<r'x-'X''xx':;:7r'x  —  rrxx,,,.,i-xr:’rxx'- — ,  r:  f  i  r  *  -r  r 

IF  (INC.LE.200)  JUMP  =  2 
I F  ( MOD  (IG-1, JUMP )  . EQ .  0 )  THEN 

KOUNT  =  KOUNT  +  1 
IF(IG.EQ.l)  THEN 
KOUNT  =  0 
GOTO  886 
END  IF 

KAE( KOUNT)  =  KAA(IG) 

886  END  IF 

IF  ( DP . GE . D ( NUM  +  1))  NUM  =  NUM  +  1 
18  CONTINUE 
C 

CALL  El GENS (BOTTOM, KOUNT, M, KAE, EIC2 , RMIM) 

C 
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*  *  -T  X  * 


•ft  JT 


EXACT 


^  r-  m 


Q:V  *  T  A'  T  t  *  ■<  X  X  X  7f  *  ^  r  r  X- 

C-  ABSORPTION  IN  NEPERS/METER  CALCULATED  AS  FUNCTION  OF  FREQ 


FREK  =  FREQ/ 1000 
FREER  =  FREK- -2 

ALFAW  =  . 00001  * FREER/ ( 1  + FREER) 

C 

Q-r.v  x.t  *  i  irx  jr  r  *  r  i  7i'  x  it  '■  ’t  x  x  jr  *  x  ft  i  -r  *  ir  *  r  x  >■  i  x  '  sScrr*)r-r7(sifr7<,A,*TT7t7<-!r*'<:'\,rx*Titx7,'x7rrTT 


C*  NEXT  FOLLOV/S  A  LOOP  WHICH  CAN  EE  USED  TO  DETERMINE  THE  PRESSURE 

C*  OR  ATTENUATION  AT  ALL  DEPTHS  FROM  TOP  TO  BOTTOM 

c 

DO  54  IDR  =  37,37 
DEPTHT  =  1.0* IDR 
DO  55  I  = 1 , M 
CRIT-CRITA 
NC  =  1 
ALFAB  =  0.0 
ALFAT  =  ALFAW  +  ALFAB 
ALFA  =  CMP LX (0.0, ALFAT) 

5500  DO  598  IK  =  1,INC*N+1 

R ( IR )  =  RAA(IK)  +  ALFA 
598  CONTINUE 

ARGUE=K(20)*-2-(PI* ( I - . 5 ) /BOTTOM ) * *2 
CALC(I)  =  CSQRT( ARGUE) 

Qx*jrif*7H:7t*'rxx7(x*7t7t*>.'+rj"<r*xxrx'<rxx**)rjtT*TxxTTTrx-ivj:*'*+Titix>i7t7i:  +  7riri'xxj:TT 

C*  CALC  IS  USED  TO  EVALUATE  THE  CORRECTNESS  OF  THE  PROGRAM 

C*  BY  INPUTTING  AN  ISOVELOCITY  PROFILE  IT  IS  FOSSIBLE  TO  CALCULATE 

C*  THE  CORRECT  HORIZONTAL  WAVENUMBER  FOR  EACH  MODE  AND  THEN  COMPARE 

C*  THE  RESULTS  OF  THE  PROGRAM  TO  THIS  'CORRECT1  ANSWER 

Q>nl;'irHr*)f7!X7nt7()r****>:*T7ri(ix,\'X7irxri*xx**7[5VX»f+)txrl,'sx**>(xS7r>,,rrtrr:7(j:x7r1'7r''  »-rxir 

VAL  =  EIG2 ( I ) 

CORRN  =  1 ■ 


IF  (TCP)  THEN 


C*  IF  THE  SOUND  VELOCITY  IS  GREATER  AT  THE  TO?  THAN  AT  THE  BOTTOM 

C*  THEN  SUBROUTINE  M0DES1  IS  USED;  OTHERWISE  MODES 2  IS  USED 

C*  AFTER  EACH  RUN  OF  EITHER  SUBROUTINE,  THE  CORRECTION  MADE  TO  THE 
C*  SQUARE  OF  THE  HORIZONTAL  WAVELENGTH  IS  COMPARED  TO  A  ERROR 

C*  CRITERION;  IF  IT  EXCEEDS  THE  SET  LIMIT  THE  SUBROUTINE  IS  CALLED 

C*  AGAIN.  IF  THE  CORRECTION  IS  LESS  THAN  THE  LIMIT,  THE  PROGRAM 

C*  PROGRESSES  AND  SOLVES  FOR  THE  PRESSURE  FOR  THAT  MODE 

Q**rA,*r*Vtx**5tit**x*^'i'**'|fiit**ilr***7ir**?filf**7r7r**i)n\**x***irf*x***J'*x*7fT***  +  *7r  +  r 

c* 

180  CALL  MODES  1 ( VAL , DELZ , M, K , P , CCRR ,  INT, INC, RATIO, CRIT, CORRN,  I 

NC  =  NC  +  1 

C  I F ( CABS ( VAL ) .GT.CABS(K( INC*N+1)**2) ) THEN 

C  CRIT  =.5*CRIT 

C  CORRN  =  CORN 

C  IF  (CRIT. LT. .005 )  GO  TO  66 

C  VAL  =  EIG2 ( I ) 

C  GOTO  180 

C  END IF 

I F ( CABS ( CORR ) . GT . ERRMAX )  GOTO  180 
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o  o  o  o  o  n  n  n 


EXACT 


7  :  IE  : 


IF  (CRIT. LT. CRITA)  THEN 
CRIT  =  CRITA 
GOTO  180 
END  IF 
GOTO  57 
WRITE ( 4 , 570 ) I 

I F  ( VAL • EQ . 0 )  THEN 
ELSE 

190  CALL  MODES  2  (  VAL,  DELE  ,  N ,  R,  ?,CGRR,  I  NT ,  INTC  ,  RATIO ,  CRIT,  CCERX.’,  I  ) 

NC  =  NC  +  1 

I  F  (  CA5S  ( VAL )  .GT.CABS(K(  INC"::+1  )  *-2  )  )  THEM 
CRIT  = . 5  *CRIT 
CORRN  =  CORN 

IF  (CRIT.LT. .005)  GO  TO  56 
VAL  =  EIG2 ( I ) 

GOTO  190 
END  IF 

I F ( CABS ( CORR ) . GT . ERRMAX )  GOTC  190 
IF  (CRIT.LT. CRITA)  THEN 
CRIT  =  CRITA 
GOTO  190 
END  IF 
GOTO  57 

56  WRITE(4, 570) I 

570  FORMAT ( '  NO  CONVERGENCE ! ! ! !  »  FIND  THE  ERROR' .16) 

5  7  ENT)  IF 

E I G ( I )  =  CSQRT( VAL) 

IDT1  =  IFIX(DEPTHT/(DELZ-2 ) )-2 

IDT2  =  I DTI  +  2 

PDIFF  =  P ( IDT 2 )  -  P(1DT1) 

DDIFF  =  DEPTHT  -  ( IDT1 ) *  DELE 

PTX  =  P(IDT1)  +  PDIFF*DDIFF/(DELZ-2 ; 

IDR1  =  IFIX(DEPTHR/(DELZ*2)  ).*2 

IDR2  =  IDR1  +  2 

PDIFF  =  P( IDR2 )  -  P( IDR1) 

DDIFF  =  DEPTHR  -  (IDR1)*DELS 

PRX  =  P(IDR1)  +  PDIFF-DDIFF.  (DELE-:) 

AIG(I)  =  CMPLX(KEAL(EIG( I )  ) , AES ( A I HAS ( EIC ( I ) ) ) ) 

DENOM  =  INT*CSQRT(AIG(I )*FI/2) 

PPRESS ( I ) =PTX~PRX/DENOH 

C  FOLLOWING  INSTRUCTIONS  ONLY  FOR  CALC  CF  RAYLEIGH  COEFFICIENT 
C  DONE  FOR  RANGE  OF  CRITICAL  ANGLES  AND  FOR  DIFFERENT  IHPEOANC 

55  CONTINUE 
C 
C 

C*  THE  FOLLOWING  LOOP  CALCULATES  THE  ATTENUATION  OF  THE  SIGMA 
C*  DIFFERENT  RANGES  BY  SUMMING  THE  COMPLEX  PRESSURE  DUE  TO  TH 
C*  INDIVIDUAL  MODES 

C*  CAN  BE  SET  UP  TO  PRODUCE  ATTENUATION  AT  REGULAR  RANGE  INTERVALS 

C*  OR  AT  IRREGULAR  INTERVALS  SUCH  AS  EVERY  10  METERS  OUT  TO  ICO 

C*  METERS , THEN  EVERY  100  METERS  TO  1000  AND  SO  ON. 
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to  M  ;  k  [i] 


o  o 


IMAG=CM?LX ( 0 , 1 ) 

DO  560  IRKG  =1,1 

DO  561  NRMG  =  300,3000 
R=(10**IRNG)>-NRNG 
C  R  =  50000 

PRESSR  =  0 
STRONG  =  0 
DO  562  I  =  1 , N 

PRESS ( I ) =P?RESS ( I ) *CEXP ( I MAG*  E I G ( I  ) -R  ) ,/SQRT(  R) 

MAG (  I  )  =  CABS ( PRESS ( I ) ) 

PRESSR=  PRESSR  +  PRESS(I) 

PR ( I ) =REAL( PRESSR) 

PIM(I)  =  A IMAG ( PRESSR ) 

EASE  =  AT AM (AI MAG (PRESSR) /REAL ( PRESSR ) ) 

STRONG=CABS ( PRESSR ) 

C  WRITE  (’4,3434)1,  PRESS  (  I  )  ,  PPRESS  (  I  ) 

3434  FORMAT ( I6,4E15.7) 

C  PRINT  3487, MAG ( I ) , STRONG 

3  487  FORMAT ( 2  F15 . 5 ) 

FACT  =  CSQRT(K(INC*N+1)**2  -  ( EIG( I )/COS (CRITA) ) **2 ) 

REF  =  (RATIO  -  FACT )  /  ( R.ATIO  +  FACT) 

FRACT  =  CABS (REE) 

COMP  =  CSQRT ( K ( INC*N+1 ) **2-EIC( I ) **2 )/(2*K( IHC*N+1 ) x  BOTTOM ) 
QLOS  =  (1 -FRACT) -CABS (COMP) 

PHAS  =  AT AN ( A I MAG ( REF ) /REAL ( REF ) ) 

THETA  =  ARCOS ( E IG ( I ) /K ( INC-N+ 1 ) ) 

C 

562  CONTINUE 

510  CONTINUE 

ATTEN  =  -  20'r  ALCG10  (  STRONG) 

561  CONTINUE 

560  CONTINUE 

54  CONTINUE 
1 1  STOP 
END 


SUBROUTINE  MODES  1  (  VAL,  DELE ,  N ,  K,  ?, CCRR,  I  N'T,  INC,  RATIO,  CR  IT,  COR.RN,  I  ) 
COMPLEX  K( 10001) , P( 10001 ) , VAL, CORE, INT, FI , A1,B1  F2 , A2 , B2 , F3 , A3 , 
1B3, F4, B, A, IMAG, AP , DENO 
C 

Cx  THIS  SUBROUTINE  USES  A  FOURTH  ORDER  RUNGS  KUTTA  TECHNIQUE 

C*  TO  CALCULATE  FRESSURE  AND  ITS  DERIVATIVE  STARTING  AT  THE  SURFACE 

C*  AND  STEPPING  DOWN  A  DISTANCE  2  TIMES  THE  DEPTH  INCREMENT  UNTIL 
C*  THE  EOTTOM  IS  REACHED 

C-  CALC  AT  EACH  L  IS  IN  FACT  FOR  THE  L+l  STEP 

C*  AT  THE  EOTTOM  THE  CALCULATED  PRESSURE,  DERIVATIVE  AND  INTEGRAL 
C-  ARE  USED  TO  DETERMINE  THE  CORRECTION  TO  BE  APPLIED  TO  THE 

C*  SQUARE  OF  THE  BEST  GUESSS  OF  THE  NORMAL  MODE  HORIZONTAL  WAVE- 
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HUMBER  THE  CORRECTION  IS  APPLIED  AND  THE  RESULT 
TO  THE  MAIN  PROGRAM 

IN  ADDITION  THE  PRESSURE  AT  ALL  DEPTH  INTERVALS  IS 
THIS  IS  USED  IN  THE  CALCULATION  OF  ATTENUATION 
SHOULD  SAVE  ONLY  THOSE  PRESSURES  REQUIRED  NOT;  ALL 


PI  =  3 . 141592654 
I NT  =  0.0 
E  =  0.0 


C 


C 


c 


c 


P(l)  =  B 

DO  18  L  =  1, INC*N  -  1 , 2 


D?  =  DELZ* ( L+ 1 ) 

FI  =  -3  *  ( K ( L ) ~ * 2  -  VAL) 

A 1  =  A  +  DELZ*F1 
Bl  =  B  +  DELZ*A 

F2  =  -Bl*  (K(  L+l )  **2  -  VAL) 

A2  =  A  +  DELZ  *  F2 
B2  =  B  +  DELZ  *A1 

F3  =  -B2*(K(L+1)**2  -  VAL) 

A3  =  A  +  2*DELZ*F3 
B3  =  3  ♦  2-DELZ*A2 

F 4  =  -B3*(K(L  +  2)**2  -  VAL) 

B  =  B  +  2*DELZ*(A  +  2*A1  +  2*A2  +  A3)/6 
P ( L+ 1  )  =  B 

A  =  A  +  2  *DELZ* ( F 1  +  2*F2  +  2*F3  +  F4)/6 
INT  =  INT  +  3**2*DELZ*2 


IS  PASS 
PASSED  ; 


v  M  + 1 


c 

18  CONTINUE 

IMAG  =  CMP LX (0.0, 1.0) 

FILL  =  REAL(K(INC*N+1) )**2 
DENC  =  CSQRT ( (FILL  -  VAL ) /FI LL ) 

AP  =  CSQRT  (  FI  LL-VAL/COS  (  CRIT  )  -■  *  2  )  I  MAG/RAT  10*  DENG 

CORR  =  B*(A+AP)/INT 
VAL  =  VAL  -  CCRR/CORRN 
21  RETURN 
END 


SUBROUTINE  MODES 2  (VAL,  DELZ  ,  N, , P , CORR,  INT,  INC, RATIO , CRIT, C 
COMPLEX  K( 10001) , P( 10001 ) , VAL, CORR, INT, FI , A1 , Bl , F2 , A2 , E2 , 

1 F3 , A3 , B3 , F4,B, A, IMAG , DENO 
C 

Q*i(iirr*fS;K***f*jM«r***x****?fx*s**i«rr  +  **x*'.\;'.Vx*)fif***K***Txx*rir*xjft**rift 

C*  THIS  SUBROUTINE  USES  A  FOURTH  ORDER  RUNGE  KUTTA  TECHNIQUE 
C*  TO  CALCULATE  PRESSURE  AND  ITS  DERIVATIVE  STARTING  AT  THE  B 

C*  AND  STEPPING  UP  A  DISTANCE  OF  2  TIMES  THE  DEPTH  INCREMENT  U 

C*  THE  SURFACE  IS  REACHED 

C*  CALC  AT  EACH  L  IS  IN  FACT  FOR  THE  L- 1  STEP 

C*  AT  THE  TOP,  THE  CALCULATED  PRESSURE,  DERIVATIVE  AND  INTEGRA 

C*  ARE  USED  TO  DETERMINE  THE  CORRECTION  TO  BE  APPLIED  TO  THE 
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C-  SQUARE  OF  THE  BEST  GUESSS  OF  THE  NORMAL  MODE  HORIZONTAL 

C-  NUMBER  THE  CORRECTION  IS  APPLIED  AND  THE  RESULT  IS  PAS 

C*  TO  THE  MAIN  PROGRAM 

Cx  IN  ADDITION  THE  PRESSURE  AT  ALL  DEPTH  INTERVALS  IS  PASSED 

C*  THIS  IS  USED  IN  THE  CALCULATION  OF  ATTENUATION 


PI  =  3 . 141592654 
I  NT  =  0 

LFLOOR  =  INC*N  +  1 
B  =  1. 

P ( IMC*N)  =  B 
I MAG  =  CMP LX (O.O, 1.0) 

FILL  =  REAL ( K ( LFLOOR ) ) *  *  2 
DENO  =  CSQRT ( (FILL  -  VAL) /FILL) 

A  =  CSQRT  (  FI LL- VAL/C05  ( CRIT )  *  *2 ) IMAG/RATI O/DENO 
A  =  0. 


IF 

(  I 

.  Cl 

o 

II 

< 

Z 

DO 

28 

L 

=  1 , INCXN- 1 , 2 

DP 

= 

DELS* ( LFLOGR-L-2 ) 

FI 

= 

-B  *  ( K ( LFLOOR- L+ 1 ) * 

x  2  -  VAL) 

A1 

= 

A  +  DEL-ZXF1 

B1 

= 

B  +  DELZ*A 

F2 

— 

-Bl* (K(LFLOOR-L)**2 

-  VAL) 

A2 

= 

A  +  DELZ-F2 

B2 

= 

B  +  DELZ-A1 

F3 

- 

-B2X ( K( LFLOOR- L) **2 

-  VAL) 

A3 

= 

A  +  2*DELZ-F3 

33 

B  -  2XDELZ*A2 

F4  =  -S3*  (K(LFL00R-L-1)XX2  -  VAL) 

B  =  B  +  2*DELZX(A  +  2*A1  +  2-A2  +  AS )/6 
? ( LFLCOR-L-2 )  =  B 

A  =  A  +  2 >'DELZf'  ( FI  -  2-F2  +  2*F3  -  F4)/6 
I NT  =  I NT  +  B**2XDELZ*2 

23  CONTINUE 

CORR  =  B*A/IMT 
VAL  =  VAL  +  CORR/CORRM 
21  RETURN 
END 


SUBROUT I NE  E I GENS ( BOTTOM , KOUNT , M , KAE , E I C2 , RM I N ) 

DIMENSION  E( 10001) ,D( 10001) , AK2 ( 10001 ) , E IG2 ( 50 ) , EIC(50) 
REAL  KAE (10001) 

C  PRINT  4356, BOTTOM, KOUNT, M, RM I N 

c* 

C*  COMPLEX  NOT  USED  IN  THIS  SUBROUTINE  BECAUSE  IT  CALLS  IMSL 
C*  SUB  EQRT1S  AMD  IT  IS  IMPRACTICAL  TO  MODIFY  THAT  FOR  COMP Li 

C*  THE  REAL  EIGENVALUES  ARE  THEREFORE  CALCULATED  AND  USED  IN 

C*  MODES 1  AND  MODES 2  WHERE  COMPLEX  EIGENVALUES  ARE  USED 
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ooooo  u  u  u  uuuu 


J  . 


UHIS  SUBROUT INF  IS  MOT 


JSED  I  Is  THE  FFT  PROGRAM 


DZ  =  EOTTOM'( MOUNT- 1) 

DZ2  =  DZ*DZ 

DO  100  1  =  1,  KOUNT 
AK2(I)  =  KAE(I)*~2  . 

D(I)  =  -AK2(I)-DZ2  +  2.0 
E(I)  =  1.0 
100  CONTINUE 

D ( KOUNT )  =  D ( MOUNT )  -  I . 0 
E( 1)  =  0.0 
ISW  =  0 
IER  =  0 

CALL  EQRT1S  (D, E, KOUNT, M, ISW, IER) 


DO  6000  J  =  1 . M 
EIC2(J)  =  -D(J)/DZ2 


ALTHOUGH  THE  FORMULA  IS  EXACT  THE  SUBROUTINE  RUNS  IN  SINGLE 
PRECISION  SO  ONE  MUST  BE  CAREFUL 

fcxivxxxxxxxxxxx'x*1,  Xxxx-J  x  :v  x-]t,\x+-TX*\,xxxx  +  xx'":-'r'A,xx*x7fxxxx'x'X'xx:i‘x'xx>ixxx:<:’x:rrr'1'x'r'x 

EIG2EX  =  AK2(2C)  -  ( 3 . I415926S35* (2* J-l )/( 2 . 0* BOTTOM) ) ' "2 
6000  CONTINUE 
C 

9000  FORMAT  (1 6,3  FI  6. 10) 

8000  RETURN 
END 
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APPENDIX  C 


COMPUTER  PROGRAM  FFP 


FILE:  FFP 


FORT 


COMPLEX  WRON,  U, DU, V, DV# VI , DV1 , K( 25000 ) , FN, EIGVAL , VAL ,  I MAG , FUNC , 
1AREA,  SIMPLE, SOLN, VEE, YEW,  VRTWN ,  WRCNK 
REAL  KMAX,D(30) , KAY (30) , RATE ( 30) ,0(30) 

C 

NUM=1 

PI  =  3 . 14159265 
CMAX=0 
CMIN=2000 
C 

FREQ  =  50 
DEPTHT  =10 
DEPTHR  =  37 
RMIN  =  40 

C  INC  ,  THE  NUMBER  OF  INCREMENTS  PER  REAL  MODE ,  MUST  BE  EVEN) 

INC  =  200 

12  READ  (3,*)  D(NUM) ,C(HUM) 

C 

IF  (D(NTJM)  .  LT.  0)  GOTO  14 
KAY(NUM)  =  2  *  P I  *  FREQ  /  C(NUM) 

IF  (C(NUM) .GT.CMAX)  THEN 
CMAX=C ( NUM ) 

DMAX=D ( NUM ) 

END  IF 

IF  ( C ( NUM ) . LT . CM I N )  THEN 
KMAX=KAY ( NUM ) 

DKMAX=D ( NUM ) 

CM IN  =  C(NUM) 

END  IF 

IF  (NUM  . EQ.  1)  GOTO  13 
C 

RATE  ( NUM- 1 )  =  (KAY (MUM)  -  KAY (MUM-1))  /  ( D ( NUM ) -D ( NUM- 1 ) ) 

C 

13  NUM  =  NUM  +  1 
GOTO  12 

14  BOTTOMED ( NUM- 1 ) 

N  =  I F I X ( 2  *  BOTTOM +  FREQ/C MAX  + . 5 ) 

INCTOT  ='  INC*N 

IF  ( INCTOT. GT. 18000)  INCTOT  =  18000 
DELZ  =  BOTTOM/I NCTOT 
IF  (N.GT.50)  THEN 
WRITE ( 4 , 100 )  N 

100  FORMAT (  ' TOO  MANY  MODES  FOR  THIS  PROGRAM!!  N ( REAL )  =  A  14) 

GOTO  11 
END  IF 


VALUE  OF  ABSORPTION  DEPENDENT  ON  FREQUENCY  ONLY-NO  TEMP  DEPENDENCE) 

FFREQ  =  FREQ/1 COO 

ALFA1  =  . 0000092 /  (  .  7  +  FFREQ**2 ) 

ALFA2  =  .0046/(6000  +  FFREQ* *2) 


nonnoon  o  o  no 


i':r 


A1 


FORTRAN 


ALEA3  =  , 000000046 

ALFA  =  ( ALFA1  ♦  ALFA2  +  ALFA3 )  *  FFREQ*  x2 
ALFA  =  .00001 
ALFA  =  0 

****  EXAGGERATE  ALFA  TO  SEE  IF  NUMBER  OF  ITERATIONS  IMPORTANT 
HUM  =  1 

DO  18  IG  =  1,I>:C*N  +  1 
DP  =  ( IG-1)  *  DELE 

K( IG)  =  KAY (HUM )  +■  RATE(NUM) * (DP-D(NUW) ) 

IF  ( DP . GE . D ( HUM  +  1))  HUM  =  NUM  +  1 
18  CONTINUE 

IF  ( DEPTHT . GT . DEPTHR )  THEN 
DUP  =  DEPTHR 
DLVJ  =  DEPTHT 
ELSE 

DLW  =  DEPTHR 
DUP  =  DEPTHT 
END  IF 

IU1  =  I  F IX ( DUP/' ( DELE *2  )  )  x2 
IU2  =  IUl  +  2 

I  LI  =  IFIX ( DLW/ ( DELE *2 )  )  *2 
IL2  =  IL1  +  2 


ALFA  IS  CALCULATED  FOR  ONE  FREQUENCY  AND  IS  PART  OF  VAL ( COMPLEX ) 

*K'Tiifrxx*r,)f*T-ifX7:'fc*X*,r-rT,*jf*xx*r7tT  +  7rTiTir,1;jrw>-?crr'x'X'/X,X7rrTXXT>:r*rxi1.,  +  r-Trrt  w  — 

AREA  =0.0 
I FFP  =  2047 
DO  55  I  =1, IFF? 

C 

c 

C  DURING  THIS  PHASE  OF  THE  PROGRAM  MODES 2  PROVIDES  VALUES  FOR 

C  PRESSURE  AND  DERIVATIVE  AT  BOTH  THE  TRANSMITTER  AND  THE  RECEIVER 

C  STARTING  WITH  THE  INITIAL  CONDITIONS  AT  THE  SURFACE 

C  (IE  PRESSURE  =  0  AND  DERI V  IS  AN  ARBITRARY  VALUE) 

C  MODES  2  PROVIDES  VALUES  FOR  PRESSURE  AND  DERIVATIVE  FOR  THE 

C  LOWER  POINT  ONLY  STARTING  WITH  THE  INITIAL  CONDITIONS  AT  THE 

C  BOTTOM  (IE  DERIVATIVE  =  0  AND  PRESSURE  IS  AM  ARBITRARY  VALUE 

C  A  VALUE  FN  REPRESENTS  THE  CONTRIBUTION  OF  EACH  ’INCREMENT' 

C  AND  THE  AREA  UNDER  THE  FN  CURVE  REPRESENTS  THE  TOTAL  ’PRESSURE 

C 

+  jrrxifjrrTr* 

c 

R  =  500 

WVNMCH  =  1 . l-KMAX/IFFP 

WVNM  =  I* WVNMCH 

EIGVAL  =  CMP LX (WVNM, ALFA) 

VAL  =  EIGVAL*  x  2 

CALL  MODES 1 ( VAL, DELE , INCTOT, K, U, DU, IUl , IU2 , DUP, DLW, VI , DV1 , 
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FORTRAN 


r I LZ :  REP 


1  IL1,IL2) 

CALL  M0DES2 ( VAL, DELZ ,  INCTOT , K, V, DV , I  LI,  IL2,DLW,  INC) 

WRON  -  V*DV1  -VI  "'DV 
878  FORMAT  (10E8.1) 

C  CALCULATE  PARTIAL  FIELD  INTEGRAL  FOR  PARTICULAR  VALUE  OF 

C  WAVENUMBER  COMPARE  TO  THAT  CALCED  BY  FRCGRAM 

+  +  xtT-rtxT'tr’"rf 

C  VRTWN  =  CSQRT(K(20)**2  -  VAL ) 

C  YEW  =  S IN ( VRTWN *D UP ) 

C  VEE  =  COS ( VRTWN* ( EOTTOM-DLW ) ) 

C  WRONK  =  VRTWM^COS ( VRTWN*  BOTTOM ) 

C  SOLN  =  YEW  *  VEE  /  WRONK 

C  SIMPLE  =  U* V/WRON 

FN  =  U* V*CSORT ( E IGVAL ) /WRON 
I MAG  =  CMP LX (0,1) 

FUNC  =  FN*CEX? ( IMAG*E IGVAL*  R )  *  SQRT ( 2/( PI +R ) ) 

CONT  =  CABS (FUNC) 

AREA  =  AREA  +  FUNC*WVNMCH 
TOTAL  =  CABS (AREA) 

IF  ( MOD (1,5). EQ . 0 )  WRITE  (4,190)  E IGVAL, CONT , TOTAL 
I F ( MOD (1,100). EQ . 0 )  PRINT  190,  E IGVAL, CONT , TOTAL 
55  CONTINUE 

C 

11  STOP 

END 

C  MODES 1  PROVIDES  PRESS  AND  DERI V  AT  BOTH  UPPER  AND  LOWER 

C  INPUTS 

C  VAL  -  THE  INCREMENTAL  HORIZONTAL  WAVENUM  -  U?  TO  8192  VALU 

C  DELZ-  DEPTH  INCREMENT  EN  OF  DEPTH,  NUMB  OF  MODES,  AND  NUMB 

C  OF  INCREMENTS  /  MODE 

C  N  -  NUMBER  OF  MODES 

C  K  -  HORIZONTAL  WAVENUMBER  (COMPLEX)  AT  EACH  DEPTH  INCREMENT 

C  DUP  -  DEPTH  OF  UPPER  OF  TXER  AMD  RXER 

C  DLW  -  DEPTH  OF  LOWER  OF  TXER  RXER 

C  IU1/IU2  INCREMENT  NUMBER  FOR  INC  ABOVE  AND  BELOW  'UPPER' 

C  IL1/IL2  INCREMENT  NUMBER  FOR  INC  ABOVE  AND  BELOW  ’LOWER’ 

C  INC  -  NUMBER  OF  INCREMENTS  OF  DEPTH 

C  OUTPUTS 

C  U1  -  ’PRESSURE’  AT  UPPER 

C  DU1  -  ’DERIVATIVE  AT  UPPER 

C  VI  -  'PRESSURE'  AT  LOWER 

C  DVl  -  ’DERIVATIVE’  AT  LOWER 

C 

c 

SUBROUT I ME  MODES 1 ( VAL , DELZ , I NCTOT , K , U , DU , I U1 , I U2 , DU? , DLW , V 1 , DV 1 , 
1IL1, IL2 ) 

COMPLEX  K( 25000) , P ( 2 ) , DP ( 2 ) , DPI ( 2 ) , P 1 ( 2 ) , FI , A1 , B1 , F2 , A2 , B2 , 

1 F3 , A3 , B3 , F4 , A , B , PD  I FF , P1DI FF , QD I FF , Q1D I FF , U , VI , DU , DVl , VAL 
C 

B  =  0.0 
A  =  .00001 
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1:1  W 


o  o  o  o  o 


FFP 

FORTRAN  A1. 

DO  18 

L 

=  1 , IMCTOT  *  1, 

- 

FI 

= 

-E  -  (K(L)** 2  - 

V  A  L ) 

A1 

= 

A  +  DELZ’’ FI 

B1 

= 

B  +  DELZ 'A 

F2 

- 

-Bl*(K(L*l)*+2 

-  VAL) 

A2 

= 

A  +  DELZ  *  F2 

32 

= 

B  +  DELZ  *  A1 

F3 

— 

-B2*(K(L+l)-^2 

-  VAL) 

A3 

= 

A  -  2~DELZ*E3 

B3 

= 

B  +  2*DELZ* A2 

F  4 

-53*(K(L+2)**2 

-  VAL) 

B  =  B  +  2*DELZ*(A  +  2*A1  +  2-A2  +  A3}/6 
A  =  A  +  2*DELZ-(F1  +  2*F2  +  2-F3  -  F4)/6 

I F ( (L  +  l) .EQ. IU1)  THEN 
P(2)  -  B 
DP ( 1 )  =  A 
END  IF 

IF  ( (L+l) .EQ. IU2)  THEN 
P  (  2  )  =  £ 

DP ( 2 )  =  A 
END  IF 

IF  ( (L+l)  .EQ.  I  LI )  THEN 
PI ( 1 )  =  B 
DP 1(1)  =  A 
END  I F 

IF  ( (L+l) .EO. IL2)  THEN 
PI (2 )  =  B 
DPI (2)  =  A 
END  IF 

18  CONTINUE 

FDIFF  =  ? ( 2 )  -  ? ( 2 ) 

P1DIFF  =  P 1 ( 2  )  -  PI ( 1 ) 

QDIFF  =  DP ( 2 )  -  DP ( 1 ) 

QiDIFF  =  DP  1(2)  -  DP 1(1) 

DDIFI  =  DUP  -  ( IU1)*DELZ 

DID I  Ft  =  DLW  -  (IL1)*DELZ 

U  =  p ( 1 )  +  PDI FF  *DDI FF/( DELL- 2 ) 

VI  =  Pl(l)  +  PlDIFF*DlDIFF/(DELZ"-2) 

DU  =  DP ( 1 )  +  QDIFF*DDIFF/(DELZ*2) 

DV1  =  DP 1(1)  +  Q1DIFF*D1DIFF/(DELZ*2) 

21  RETURN 
END 


•fr'i'-K-k-k-k'le-X’k-k-kTr-x-kTr-k-x'k-k  rr 

MODES2  PROVIDES  PRESS  AKD  DER  AT  'LOV.'ER'  ONLY 
INPUTS 

VAL  -  THE  INCREMENTAL  HORIZONTAL  WAVEN'JH  UP  TO  819 
DELZ-  DEPTH  INCREMENT  FM  OF  DEPTH,  NUMB  OF  MODES, 


94 


FILE;  EtrP  FORTE.*::  El 


C  OF  INCREMENTS  /  MODE 

C  N  -  NUMBER  CF  MODES 

C  K  -  HORIZONTAL  WAVENUMBER  (COMPLEX)  AT  EACH  DEPTH  INC 

C  DUP  -  DEPTH  OF  UPPER  OF  TXER  AND  RXER 

C  DLW  -  DEPTH  OF  LOWER  OF  TXER  RXER 

C  IU1/IU2  INCREMENT  NUMBER  FOR  INC  ABOVE  AND  BELOW  1 U?P 

C  IL1/IL2  INCREMENT  NUMBER  FOR  INC  ABOVE  AMD  BELOW  'LOW 

C  INC  -  NUMBER  OF  INCREMENTS  OF  DEPTH 

C  OUTPUTS 

C  U1  -  'PRESSURE'  AT  UPPER 

C  DU1  -  'DERIVATIVE  AT  UPPER 

C  VI  -  'PRESSURE'  AT  LOWER 

C  DV1  -  'DERIVATIVE'  AT  LOWER 

C 

c 

SUBROUTINE  MODES2 ( VAL , DELZ ,  INCTOT, K, V, DV,  I  LI,  IL2,DLW) 
COMPLEX  K ( 25000 ) , P ( 2 ) , DP (2)/Fl/Al,51,F2,A2,B2, 

1F3  ,  A3 , B3 , F4 , A , B , PDI FF , QD I FF , V , DV , VAL 
B  =  .00001 
A  =  0.0 


DO  28 

L 

=  1 ,  INCTOT-1- 1  ,  2 

LI 

= 

INCTOT  -  L  -  1 

FI 

= 

-B  *  (K(Ll)-*2 

-  VAL) 

A1 

= 

A  +  DELZ*  FI 

Bl 

= 

B  +  DELZ-A 

F2 

= 

-Bl*(K(Ll-l)*-2 

-  VAL) 

A2 

= 

A  +  DELZ*  F2 

B2 

= 

B  +  DELZ*Ai 

F3 

- 

-B2*(K(L1-1 )**2 

-  VAL) 

A3 

= 

A  +  2  *DELZ*  F3 

B3 

= 

B  +  2  *DELZ* A2 

F-l 

-B3 * ( K ( LI -2 ) *  *  2 

-  VAL) 

B  =  B  +  2’-DELZ-(A  *  2*A1  -  2-A2  +  A3)  '6 
A  =  A  +  2*DELZ~  ( Fi  +  2*-F2  +  2^F3  *  F4)/6 
C 

I F ( (Ll-2)  .E Q. I  LI )  THEM 
P(l)  =  B 
DP(1)=  -A 
END  IF 

I F ( ( LI -2  )  . EQ .  I L2 )  THEN 
P(2)  =  B 
DP ( 2 )  =  -A 
END  IF 
C 

28  CONTINUE 

PD IFF  =  P ( 2 )  -  ?(1) 

QD IFF  =  DP ( 2 )  -  D? ( 1 ) 

DDIFF  =  DLW  -  (ILl)-DELZ 
V  =  P(l)  +  PDIFF*DDIFF/(DELZ*2 ) 

DV  =  DP ( 1 )  +  QD I FF*DDI FF/ ( DELZ*  2 ) 

C 

21  RETURN 
END 
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